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PREFACE 

Two  feedback  control  schemes  which  maximize  a  terminal  quantity  while 
satisfying  specified  terminal  conditions  are  discussed  and  compared.   The 
schemes  are  "based  on  a  linear  perturbation  from  a  nominal  non-optimal  path 
which  does  not,  in  general,  satisfy  the  terminal  conditions.  The  methods 
have  been  programmed  in  the  ALGOL  computer  language  for  evaluation  and  the 
programs  are  included  in  the  appendices. 
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I.   INTRODUCTION 

In  recent  years  several  authors  have  treated  the  problem  of  deter- 
mining optimum  control  programs  for  nonlinear  systems  with  terminal 
constraints.  These  problems  arise  in  the  design  of  control  systems  and 
development  of  guidance  laws  where  it  is  desired  to  determine,  out  of  all 
possible  time  histories  of  the  control  variables,  the  one  control  history 
that  maximizes  (or  minimizes)  one  terminal  quantity  or  cost  function  while 
simultaneously  yielding  specified  values  of  certain  other  terminal 
quantities . 

The  steepest-ascent  or  gradient  method  developed  by  Kelley  ,  Bryson 

2 
and  Denham  ,  which  is  a  systematic  and  rapid  numerical  procedure,  has 

proved  to  be  successful  in  solving  this  class  of  problems.   Improvement 

in  the  convergence  time  of  the  iterative  process  involved  has  been  achieved 

3  k 

by  Rosenbaum  by  a  method  based  on  the  earlier  work  of  Bryson  and  Denham  . 

Another  successful  method  developed  by  Breakwell,  et. al.  ,  and  modified 
by  Bullock  is  a  second  variation  method  in  the  Calculus  of  Variations. 

The  principal  objectives  of  this  thesis  are  to  develop  a  simpler 
steepest-ascent  program  which  will  be  understandable  to  the  control 
engineer  without  a  background  in  the  Calculus  of  Variations  and  to  compare 
the  results  and  speed  of  convergence  with  the  method  developed  by  Bullock. 
In  this  way  it  is  expected  that  the  steepest-ascent  program  will  prove  to 
be  a  useful  instrument  in  education  and  research,  while  at  the  same  time 
through  the  comparison,  illustrate  the  advantages  and  disadvantages  of  the 
two  approaches  to  the  problem. 

The  method  used  in  developing  the  steepest-ascent  program  is 
essentially  a  variation  of  the  methods  of  Bryson  and  Denham  ,  and 
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Bosenbaum  ,  hence  it  is  restricted  to  problems  in  which  only  the  deviations 
in  the  control  variables  and  adjustable  parameters  are  considered  in  the 
performance  index.   It  is  also  restricted  to  problems  in  which  the  pay-off 
quantity  is  a  function  of  the  terminal  value  of  the  states.  This  varia- 
tion includes  a  terminal- error  control  scheme  which  maintains  a  bound  on 
the  terminal  constraint  errors,  hence  reducing  the  total  number  of  itera- 
tions required  to  converge  to  the  optimum  since  larger  deviations  from  the 
nominal  trajectory  can  be  tolerated  while  still  meeting  the  desired  terminal 
conditions . 

A  numerical  example  is  given  of  a  rocket  ascent  trajectory  into  a 
circular  orbit  of  maximum  altitude.  Provision  is  made  for  a  two-stage 
rocket  with  optimization  of  the  inter-stage  coast  duration. 


II.   STATEMENT  OF  THE  PROBLEM 


Given  a  system  which  can  be  described  by  a  set  of  non-linear  (or 
linear)  first  order,  ordinary  differential  equations,  determine  a  control 


history  u(t),  in  the  interval  t„  ^  t  <   T,  to  maximize 


t>  =   x1(T), 


(1) 


subject  to  constraints 


\|/  =   ^[x(T)]    =  0, 


ax 

dt 


=  f[x(t),u(t),t], 


tQ,   T,    and     x(t   )      given: 


(2) 
(3) 


where 


u 


(t)   = 


\M 


u    (t) 

m 


an  m  x  1  matrix  of 
'   control  variables, 


(5) 


x(t)  = 


xx(t) 


X  (t) 

n 


an  n  x  1  matrix  of  state 
'  variables, 


(6) 


a  q  x  1  matrix  of  terminal 
constraint  functions,  each 
J   of  which  is  a  known  function 
of  x(T), 


an  n  x  1  matrix  of  known 
functions  of  x(t),  u(t), 
and  t,  i.e.,  the  system 
equations; 


(7) 


(8) 


f)     is  the  pay-off  quantity  and  is  one  of  the 

states,  namely  x  (t);  (9) 

T  is  the  terminal  value  of  the  independent  variable. 

On  each  iteration  it  is  desired  to  minimize  the  mean  value  of  a 
positive  definite  quadratic  form  in  the  control  variable  deviations: 

T 

C  =  |  f    5u'(T)  6u(T)   dT  (10) 

2  ^t0 

where  the  superscript  (  ')  indicates  the  transpose  of  a  matrix,  and  Su(t) 
are  small  deviations  in  the  control  history  from  a  nominal  non-optimum 
trajectory. 


III.   STEEPEST -ASCENT  METHOD  OF  SOLUTION 

A .  BACKGROUND 

i 

Bryson  and  Denham,  in  Ref .  p ,  considered  terminal  control  of  non-linear 

(and  linear)  systems  for  minimum  mean  values  of  a  positive  definite  quad- 
ratic form  in  the  control  variable  deviations.  That  is,  it  was  assumed 
that  a  nominal  control  history  had  been  determined  which  caused  the  vehicle 
to  arrive  at  the  terminal  point  with  desired  values  of  certain  specified 
terminal  conditions.   Small  deviations  from  this  nominal  trajectory  were 
considered  which  might  he  caused  by  disturbances,  inaccuracies  in  the  ■ 
data,  inaccuracies  in  the  control  system,  etc.  The  problem  was  to  determine 
small  deviations  from  the  nominal  control  so  that  the  terminal  constraints 
would  be  satisfied  in  spite  of  the  disturbances. 

In  the  present  paper  the  nominal  trajectory  is  determined  by  guessing 
a  reasonable  control  variable  program.  For  example,  in  a  rocket  trajectory 
problem  one  might  choose  an  initial  launch  angle  and  a  gravity  turn  with 
zero  thrust  angle  throughout  as  is  done  in  the  numerical  example.   Further- 
more, it  is  desired  not  only  to  determine  control  deviations  which  result 
in  meeting  the  terminal  constraints,  but  also  to  maximize  the  terminal 
value  of  one  of  the  states  while  minimizing  a  performance  index.  This 
optimization  scheme  is  a  variation  of  the  so-called  Lambda  Matrix  Control 
feedback  method  described  in  Ref.  /  and  the  convergence  method  of  Ref.  %• 

B.  DERIVATION  OF  EQUATIONS 

The  optimum  programming  problem  can  be  solved  systematically  and 
rapidly  on  a  high  speed  digital  computer  using  the  steepest-ascent 
technique.  As  stated  in  (A) ,  this  technique -starts  by  guessing  a  nominal 
control  variable  program,  u*(t),  and  solving  the  set  of  differential 
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equations  (3)  with  initial  conditions  {k) ,   to  determine  a  nominal  tra- 
jectory. This  trajectory,  in  general,  will  neither  maximize  f>   nor  will  it 
satisfy  the  terminal  constraints  (2). 

Consider  small  perturbations  in  the  control  variables,  6u(t),  about 
the  nominal,  where 

6u(t)  =  u(t)  -  u*(t)  (11) 

(The  superscript  (*)  indicates  terms  evaluated  along  the  nominal  trajectory.) 
These  perturbations  will  cause  perturbations  in  the  state  variables,  5x(t), 
where 

6x(t)  =  x(t)  -  x*(t)  (12) 

Substituting  these  relations  into  the  system  differential  equations  (3)  and 
expanding  f  in  a  Taylor  series  about  the  nominal  the  result  is,  to  first 
order 

~  (5x)  =  F(t)6x(t)  +  D(t)8u(t)  (13) 

where 

an  n . x  n  matrix  of 
y(f\   _  C^\  partial  derivatives 

■  ~  \(§x/  '   evaluated  along  the 
nominal  trajectory; 

and 


D 


w  ■  (ir 


> 


an  n  x  m  matrix  of 
partial  derivatives 
evaluated  along  the 
nominal  trajectory 


To  determine  the  effects  upon  the  terminal  conditions  f)   and  ty   we 
introduce  the  linear  differential  equations  adjoint  to  (12)  defined  as 

§  (T,t)  =  -*(T,t)F(t)  {Ik) 


where  0  is  an  n  x  n  fundamental  or  state  transistion  matrix  whose  elements 
give  the  sensitivities  of  the  terminal  states  to  perturbations,  Bx(t), 
along  the  trajectory.   (See  Ref.  7) •   Initial  conditions  for  these  equations 
are  specified  at  the  terminal  time,  i.e., 


o(T,l)  =  I,   the  identity  matrix 


(15) 


hence  numerical  integration  proceeds  backward  in  time. 

The  solution  to  (1*0  provides  a  solution  to  the  linear  perturbation 
equations  (12)  at  the  terminal  point: 


5x(T)  =  <D(T,t)dx(t)  +  f       0(T,r)D(T)5u(T)dT 


(16) 


Since  ft   and  \|/  depend  on  the  terminal  values  of  the  states,  small 
deviations,  af)   and  d4  may  be  calculated  from  the  solution  to  (l6).   Con- 
versely, if  5x(t)  is  known  by  specifying  values  of  dj#  and  d\|r,  the  corres- 
ponding control  history  may  be  calculated,  which  is  what  is  done.  Re- 
writing (l6) 


Sx 


(T)  -  $(T,t)Bx(t)  -  f      $(T,r)D(T)5u(T)dT 


(17) 


In  order  to  minimize  the  performance  index  subject  to  the  constraints 
(l6),  the  method  of  Lagrange  multipliers  is  employed.  Multiplying  (l6)  by 
a  matrix  of  Lagrange  multipliers 


v  = 


LVi. 


(18) 


(where  (l6)  is  written  to  include  only  those  states  which  appear  in  'fi   and 
\[<,  hence  q+1  equations)  and  adjoining  the  result  to  (10) 

C  =  v'Sx(T)  -  v*$(T,t)5x(t)  +J        [I  5u'(t)6u(t)  -  v  ,0(T,t)d(t)5u(t)1  dr 

t  (19) 

Completing  the  square 

C  =  v'[Bx(T)  -  0(T,t)5x(t)]  +1     ||  (Su(t)D'(t)<J>'(T,t)v)  f  (6u(t)-D!  (t)$ '(T,t)v) 

-  |  v,0(T3T)D(T)D,(T)$'(T,T)vldT  (20) 

To  minimize  C  subject  to  the  control  variations  choose 

6u(t)  =  D,(t)$'(T,t)v  (21) 

Substituting  this  relation  in  (17)  in  order  to  find  v 


5x(T)  -  0(T,t)Bx(t)  -  /    $(T,T)D(T)D'(T)$(T,T)vdT  =  0  (22) 


Define  the  controllability  matrix 

rT 

J  =  +/   $(T,T)D(T)D'(T)0T(T,T)dT  (23) 

Note:   This  integration  may  be  performed  simultaneously  with  (lh)   by 
numerically  solving 

|U  -»(T,t)D(t)D»(t)*'(T,t)  (24) 

with  initial  conditions 

J(T)  =  0  (25) 

Equation  (22)  may  be  written 


6x(T)  -  <2>(T,t)8x(t)   =  Jv  (26) 

Solving  for  v 

v   =   J-1[5x(T)    -   <D(T,t)6x(t)]  (27) 

Substituting  this  relation  in   (21) 

6u(t)   =  D,(t)0,(T,t)j"1[8x(T)    -   $(T,t)6x(t)]  (28) 

which  is  the  perturbation  in  the  control  history  which  satisfies  the  con- 
straint (17)  while  minimizing  the  performance  index  (10). 

As  stated  earlier,  8x(t)  is  determined  from  the  specified  values  of 
dcp  and  d\|/  where 

dcp  =  Bac1(T)  (29) 

and 

d;  =  \l»[8x(T)]  (30) 

As  yet,   nothing  has  been   said  as  to  how  one   chooses  the  desired  pay- 
off improvement,    dcp,   or  the   desired  improvement   in  meeting  the  terminal 
constraints,   d\|/.      The   latter   is  normally  chosen  as 

d\|r  =    -ty  (31) 

that  is,  the  negative  of  the  total  error  on  any  iteration  is  chosen  as 
the  desired  correction  specified  on  the  following  iteration.  The  problem 
of  specifying  dcp  is  more  complex  and  is  the  subject  of  the  next  section. 

It  is  worthy  to  note  at  this  point  that  the  Lagrange  multipliers, 
which  are  error  feedback  terms >  need  not  be  computed  at  every  point  on 
the  trajectory.   Sufficient  accuracy  can  be  obtained  in  computing  the 
control  deviations  (21)  by  solving  (27)  at  discrete  intervals  and  using 
the  result  until  the  next  "sampling  time".   This  reduces  the  number  of 
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times  the  controllability  matrix  must  be  inverted  per  iteration  and  materi- 
ally improves  the  running  time  of  the  program.  Experimentation  will  re- 
veal how  large  the  sampling  interval  can  be  made.   Since  the  controllability 
matrix  is  singular  at  the  terminal  time,  T,  new  values  of  the  Lagrange 
multipliers  should  not  be  computed  too  close  to  the  end  of  the  trajectory. 

C.   METHOD  OF  SPECIFYING  THE  IMPROVEMENT  IN  PAY-OFF 

In  general,  one  does  not  know  how  far  from  the  optimum  a  given  nominal 
trajectory  will  be.   It  is,  therefore,  difficult  to  guess  how  much  pay-off 
improvement  to  specify  initially.  However,  it  is  possible  to  compute  a 
value  of  dcp  which  will  result  in  a  trajectory  that  satisfies  the  terminal 
constraints.   This  is  done  as  follows: 

The  changes  in  the  control  variables  required  to  meet  the  terminal 
constraints  with  the  pay-off  unconstrained  are  given  by 


Bu(t)  -  A 

where 


I        A  A'dT 


1-1 

d\|r  (32) 


A  =  $(T,t)D(t),  without  row  1  or  column  1. 

This  is  the  same  as  the  basic  control  equation  (28)  with  px(tn)  equal  to 
zero  and  5x(t)  containing  only  the  terminal  constraint  terms. 

The  change  in  pay-off,  dcp,  that  will  be  produced  by  a  given  change 
in  the  control  variables  is,  from  (l6) 


/T 
A1Su(T)dT  (33) 

-0 


whe 


re  A,  is  the  first  row  of  the  A  matrix.   Substituting  (33)  in  (32) 
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where  A-.    is  the  first  row  of  the  A  matrix.      Substituting   (33)    in  (32) 

-1 


dcp  =   6x   (T)   = 


fT  A     A'dT        fT  A  A'dr 

J*0  Jt0 


d^  (34) 


Equation  (3^)  gives  the  change  in  pay-off  associated  with  adjusting 
the  control  in  order  to  meet  the  terminal  constraints.  This  value  is  used 
on  the  first  iteration. 

Equation  (5k)   is  also  used  to  compute  a  value  of  the  pay-off  cor- 
rected for  terminal  errors,  i.e.,  cp  +  dcp  is  the  value  the  pay-off  would 
have  achieved  had  the  terminal  errors  "been  zero.  Hence,  one  can  determine 
whether  an  improvement  in  the  pay-off  was  actually  achieved  or  if  an 
apparent  improvement  was  a  result  of  larger  terminal  constraint  errors. 

On  subsequent  iterations,  one  of  three  methods  is  used  to  compute  dcp. 
A  value  equal  to  25  per  cent  of  the  nominal  value  of  cp  is  computed  and 
stored.   This  quantity  is  called  dcp**  and  is  used  in  method  (2)  below.   It 
is  a  fairly  arbitrary  choice  but  should  be  made  as  large  as  seems  reason- 
able. The  program  will  automatically  adjust  it  if  it  is  too  large. 

Method  1.   Choose  dcp  to  satisfy  the  terminal  constraints 

with  the  pay-off  unconstrained  as  described  above. 

Method  2.   If  I  d\j/ 1  ^  e,  where  s  is  chosen  as  reasonable  tolerance  on 
the  terminal  constraints  then 

a*  -  SSSL  (35 ) 

21 

where  i  is  a  count  of  the  number  of  times  method 
(3)  has  failed.   This  has  the  effect  of  halving 
the  improvement  specified  each  time  a  run  is 
unsuccessful  in  improving  the  terminal  errors  or 
the  corrected  value  of  the  pay-off.  The  program 
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terminates  when,  while  executing  this  method, 
dcp  becomes  less  than  a  pre-set  number.  A 
final  run  is  then  made  using  method  (l). 

Method  3.   If  |  di|/|  >  e  the  following  questions  are  asked: 
Were  the  errors  on  the  current  iteration 
smaller  than  on  the  preceeding  iteration? 
Was  there  an  improvement  in  the  corrected  value 
of  the  pay-off?   If  the  answer  to  either 
question  is  no,  the  run  is  considered  unsuc- 
cessful, the  control  history  is  replaced  by  the 
previous  control  history,  and  method  (l)  is 
used.   If  the  answer  to  either  question  is  yes, 
dcp  is  set  equal  to  zero  and  an  attempt  is  made 
to  satisfy  |dtl  -  £•   If  this  test  fails  a 
second  time,  method  (l)  is  used. 

D.   TWO-STAGE  ROCKET  TRAJECTORY  WITH  COAST  PERIOD  OPTIMIZATION 

In  raany  orbit  injection  applications,  such  as  the  Gemini-Titan  II 
system,  the  launch  vehicle  is  made  up  of  two  powered  stages.   It  is  there- 
fore of  interest  to  consider  the  effect  of  an  interstage  coast  phase  on 
the  maximum  altitude  obtainable.   In  this  section  a  method  of  calculating 
the  optimum  coast  duration  is  derived. 

The  basic  equations,  (l)  through  (12),  are  the  same.  The  linear 
perturbation  equations  (13)  may  be  written 

|^  [5x(t)]  -  F(t)5(t)  +  D(t)5u(t)  +  B(t)8c  (36) 


where 


an  n  x  1  matrix  of  partial 
B(t)  =  (|£)*   ,  derivatives  of  f  with  (   } 

*  '        \  oc  /     respect  to  the  coast 
duration,  c. 
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The  solution  to  (37)  is 


/T  y-rp 

0(T,T)D(T)5u(T)dT  +/    $(T,T)B(T)6cdT 


(38) 


The  performance  index  becomes 


C  =  |  /    6u'(T)5u(r)dT  +  -  A  5c         .  (39) 

where  A  is  simply  a  weighting  factor. 

Before  attempting  to  minimize  (39)  subject  to  the  constraint  (38), 
a  method  must  "be  derived  to  evaluate  the  last  term  in  (38)  which  is  the 
change  in  the  terminal  values  of  the  states  due  to  a  change  in  the  coast 
duration.   Since  c  is  an  adjustable  parameter  which  does  not  appear  ex- 
plicitly in  f ,  the  partial  derivatives  cannot  be  evaluated  directly. 
However  the  desired  term  may  be  calculated  by  the  following  method: 
Define 

t,  =  Stage  I  burnout  time 

t  =  Stage  II  ignition  time 
hence 

c  =  t2  -  tx  (40) 

Since  5c  is  a  small  time  increment  we  may  write 

t2+At 
x(t2  +  At)  -  x(t2)  =  J  g  dT 

t2 

>:(t2)+Ax 
"  x(t2) 

=  (x(t2)  +  Ax)  -  x(t2); 
13 


consequently  we  may  make  the  approximation 

x(ts  +  6=)  -  x(t2)  2  (||)t2  so  (41) 

Now  define  x  as  the  states  evaluated  with  the  thrust  off  (uncontrolled), 

u  ' 

and  x  as  the  states  evaluated  with  the  thrust  on  (controlled).   From  (hi) 
xu(t2  +  6c)  =  x(t2)  +  xu(t2)5c 
xc(t2  +  6c)  =  x(t2)  +  xc(t2)6c 

where  (')  indicates  differentiation  with  respect  to  time.   Subtracting  the 
above  expressions  we  have 

xu(t2  +  6c)  -  xc(t2  +  6c)  =  [xu(t2)  -  xc(t2)]6c 


Define 


6xc  =  [W  "  W]5C   '  (1+2) 


The  quantity  5x  is  a  perturbation  in  the  states  occurring  at  time 

t  =  t  +  6c 

due  to  cutting  off  the  thrust  for  a  period  6c.  The  question  remains: 
How  does  this  perturbation  propagate  to  the  end  of  the  trajectory?   The 
answer  is  clearly 

5x(T)  =  $(T,t  +  6c)6xc  (43) 

Finally,  (38)  may  be  written 

6x(T)  =  <5»(T,t)6x(t)  +  /   ^(T,T)D(T)6u(T)dT  +  $(T,t  +  6c)6x  (kh) 


where  the  last  term  is  zero  prior  to  time  tp  *+  6c. 


Ik 


Introducing  the  Lagrange  multipliers,   v,   and  adjoining   (44)   to   (39) 


/T 
8u'(T)6u(T)dT  +  |  A5c2  +  v*6x(T)  -  5  f$(T,t)6x(t)  (1+5) 


J         0(T,T)D(T)6u(r)dT  -  v  '0(T,t2  +  5c)6x 


'  t 
As  before,  completing  the  square  yields  the  optimum  control  change 

6u(t)  =  D'(t)0'(T,t)v  (46) 

By  differentiation  the  optimum  coast  change  is 

5c  =|v'  »(T,t2  +  8c)[xu(t2)  -  xc(t2)]  (47) 

Substituting  (46)  and  (47)  into  (44) 

•T 


5x(T)  =  Q(T,t)6x(t)  +  /   0(T,T)D(T)D,(.T)$'(T,T)dTv 

+  i  v'$(T3t2  +  6c)[xu(t2)  -  xc(t2)][xu(t2)-xc(t2)]'  $'(T,t2+Sc)v 

Define 

A  =  5x(T)  -  $(T,t)5x(t)  (49) 

J  =  f      $(T,T)D(r)D'(r)c(T,T)dT  (50) 

A  =  [xu(t2)  -  xc(t2)]  (51) 

G  =  $(T,t2  +  5c)A  A'  »'(T,t2  +  8c)  (52) 

RevTriting   (48) 

A  =   Jv  +  j-  Gv  (53) 

v  =   (j  +  f  )_1A  (54) 
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Finally,    substituting  in   (h6)   and   (V7),   the   control  and  coast  variations  are 
5u(t)    =  D'(t)C'(T,$)    (j  +  |)~     A  (55) 

5c   =  ~  A*<£(T,t2+5c)(j  +  |)     [8x(T)    -   $(T,t2  +  6c)6x(t2  +  8c)  ]  (56) 


Since  the  last  term  in  (hk)    is  zero  prior  to  time  t  +  5  ,  the  term 
G/A  in  (55)  and  (56)  is  also  zero  prior  to  that  time. 

Equation  (56)  should  he  evaluated  at  t  =  tp  +  5c,  but  since  5c  is  the 
unknown,  it  is  evaluated  at  t  .   This  does  not  introduce  an  appreciable 
error  if  6c  is  small. 

In  some  numerical  integration  procedures  the  terminal  value  of  the 
independent  variable  cannot  be  changed  once  the  integration  has  begun, 
hence  the  change  in  coast  time  cannot  be  added  immediately.  This  problem 
is  solved  by  evaluating  (56)  on  each  forward  integration  (except  the  nominal) 
and,  if  6c  f   0,  re-integrating  the  latter  portion  of  the  trajectory  from 
t,  to  T. 

E.   COMPUTATIONAL  PROCEDURES 

As  stated  earlier,  this  steepest-ascent  technique  requires  the  use  of 
a  high  speed  digital  computer.   The  sequence  of  operations  is  summarized 
here. 

1.  Compute  the  nominal  path  by  integrating  the  system  differential 
equations  (3)  with  a  nominal  control  history  and  appropriate  initial  con- 
ditions and  store  the  time  history  of  the  state  variables  at  reasonably 
snail  intervals.  Print  out  the  values  of  cp  and  1]/. 

2.  Integrate  the  adjoint  differential  equations  (l^-)  backward, 
evaluating  the  partial  derivatives  on  the  nominal  path  by  reference  to 

the  states  stored  in  step  (l).   Simultaneously  integrate  the  controllability 
matrix  equations  (2^+).   Store  the  results  at  the  same  interval  as  the  states. 
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3.   Select  desired  terminal  condition  changes,  dcp  and  d\|/,  as  ex- 
plained in  Sections  B  and  C. 

k.      Compute  and  use  the  new  control  history  while  integrating  the 
system  differential  equations  forward.  Again  print  out  the  values  of  cp 
and  v?  unless  the  next  step  applies. 

5«   If  the  two-stage  rocket  problem  is  being  solved,  compute  the  new 
coast  period  in  step  (^).  Transfer  the  storage  locations  of  the  second 
stage  control  history  to  correspond  with  the  new  coast  time.   Integrate  the 
system  equations  from  t,  to  the  new  terminal  time.  Print  out  the  values  of 
cp  and  ty. 

6.  Repeat  procedures  (2)  through  (5)  until  the  pay-off  improvement 
in  step  (3)  is  less  than  a  preset  value.  At  this  point,  use  method  (l) 
described  in  Section  C  to  select  dcp  and  complete  step  (h)   and  (5)«  This 
has  the  effect  of  eliminating  any  remaining  errors  in  the  terminal 
constraints . 

7.  Punch  cards  or  store  the  control  history  on  tape  and  terminate 
the  program. 

Before  concluding  this  section,  a  few  general  factors  of  great 
importance  in  this  type  of  numerical  calculation  should  be  discussed. 

The  programmer  must  exercise  great  care  when  working  with  values  of 
type  real  (or  floating  point).   Often  a  calculation  is  made  where  the 
result  is  expected  to  be  an  integral  value  such  as  \J2   =  2.000  ..., 
however,  due  to  the  binary,  octal,  and  decimal  conversions  which  take 
place  within  the  computer,  the  result  may  come  out  1.9999* ••99*  This 
problem  occurs  when  trying  to  generate  array  storage  indices  based  on  a 
value  of  the  independent  variable  which  is  a  floating  point  quantity. 
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IV.   SECOND  VARIATION  METHOD  OF  SOLUTION 

A.   OUTLINE  OF  THE  METHOD 

Bullock  has  derived  a  feedback  control  scheme  "based  on  the  second 
order  variational  theory  in  the  Calculus  of  Variations.  The  method  is 
outlined  here  in  sufficient  detail  to  solve  the  problem  stated  in  II. 
The  differential  equations  to  be  satisfied  are 

x  =  f(x,  u,  t)  (57) 

*  ■  "  (I)'  *  (») 


M\   /f   -Q\/M 

,N/    \-S   -F7  \N, 


and 


(59) 


b  =  M'd  -  N'c  (60) 

in  order  to  maximize 

<P  =  <P[x(T),  T]  (61) 

and  satisfy  the  terminal  constraints 

i|/  =  ttx(T),  T]  (62) 

Define  the  variational  Hamiltonian 

H  =  A'f  (63) 

The  elements  of  Eqs .  (59)  and  (60)  are 

F  =  f   -  f  H-1  H     .  (6k) 

X    u  uu  ux 

Q  =  f  H"1  f '  (65) 

u  uu 
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S  =  H   -  H   H"1  H*  (66) 

XX    xu  uu  u 


c  -  -  f  H_1  H'  (67) 

u  uu  u 

d  =  H   H"1  H'  (68) 

xu  uu  u- 

where  the  subscripts  indicate  partial  differentiation  in  the  usual  sense. 

The  initial  conditions  for  (57)  are  given  and  the  terminal  conditions 
for  (58),  (59),  and  (60)  are 

A(T)  =  (<p  -  V'*)  (69) 

x      x  t=T 

where  the  components  of  the  column  vector  v  are  sensitivities  of  the  pay- 
off <p  to  changes  in  the  terminal  constraints  ^; 

M(T)  =  I  -  Wi|r  V)"1  *    •  (70) 

.X.    J\.         J£.  J\. 

where  I  is   the   identity  matrix; 

N(T)   =   -t^tx  (71) 

h(T)    =   -tx  ^  (72) 

X'/here  d\[r  =   -ty . 

The  perturbations  in  the  control  history  are  given  by 

6u(t)  =  u*  -  u  =  -H_1(H  +  H  5x  +  f «  6A)  (73) 

v  '  uuv  u    ux      u  - 

where 

5A  =  (M')_1(N'5x  +  b).  (74) 

In  order  to  minimize  a  performance  index 
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C  =  |  /   5u'(t)5u(t)cLt  (75) 

H   in  the  above  equations  is  modified  by  adding  an  arbitrarily  large 
uu 

negative  constant,  K,  which  is  reduced  in  magnitude  as  the  program  con- 
verges.  This  has  the  effect  of  constraining  the  magnitude  of  6u. 

Since  the  terminal  conditions  on  the  adjoint  equations,  A(t),  depend 
initially  on  the  choice  of  V,  a  method  is  given  which  will  improve  the 
accuracy  of  these  terminal  conditions  on  subsequent  iterations. 

Equation  (7^0  is  an  expression  for  SA  at  any  time,  t,  but  since  M  is 
singular  at  T,  it  cannot  be  evaluated  directly.  However,  if  a  point  (t, ) 
is  chosen  sufficiently  far  from  this  singularity,  the  following  equation 
can  be  integrated  from  t-,  to  T: 

5A*  -  -S6x  -  F*  5A  (76) 

with  initial  conditions 

5A(tJ  =  (M')"1^'  Sx  +  b)  (77) 

1  t=ti 

The  solution  to  (76)  is  then  added  to  the  current  values  of  A(t)  prior  to 
the  next  backward  integration  of  (58)' 

B.   COMPUTATIONAL  PROCEDURES 

As  in  the  steepest-ascent  method,  this  method  requires  the  use  of  a" 
digital  computer.  The  sequence  of  operations  is  summarized  here. 

1.   Select  a  nominal  control  history  and  initial  values  of  v.  This 
can  be  done  by  starting  with  a  control  program  and  V  generated  by  the 
Steepest-ascent  method,  where 
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V1 

v  -  (  /   A,  A  '  dT  W  /   A  A  •  dT  ] 


(/V'*)(/ 


t 

from  Eq.  (35). 

2.  Integrate  the  system  differential  equations  as  in  the  steepest- 
ascent  method. 

3.  Integrate  Eqs.  (58),  (59)?  and  (60),  backward  with  the  appropriate 

terminal  conditions.   If  the  determinant  of  M  changes  sign  or  H   "becomes 

to     °      uu 

positive,  store  the  current  value  of  the  time  and  stop  the  integration. 
The  reason  for  this  is  explained  "below. 

k.     From  Eq.  (73),  compute  the  new  control  history  while  integrating 
the  system  equations  forward. 

5.  Compare  the  value  of  cp  and  \Jr  obtained  to  those  obtained  on  the 
nominal  trajectory.   If  the  pay-off,  cp,  or  the  terminal  constraints,  \|/, 
have  become  worse  the  run  is  considered  unsuccessful  and  a  tighter  bound  is 
placed  on  5u  by  increasing  the  magnitude  of  K.   If ,  on  the  other  hand,  cp 
and  y  are  the  same  or  have  improved,  the  run  was  successful  and  (3)  and  (h) 
are  repeated. 

6.  The  program  is  terminated  when  no  change  occurs  in  the  pay-off 

or  the  constraints  and  Jk|  «  |H   | .  At  this  point  the  control  history  is 
stored  on  punched  cards  or  tape. 

7.  If  in  step  (3)  the  determinant  of  M  changed  sign  (this  condition 
is  called  a  "conjugate  point"  in  the  Calculus  of  Variations),  or  H 
became  positive  (which  indicates  the  Legengre  condition  is  not  satisfied), 
the  integration  in  (h)   is  begun  at  a  slightly  later  time  than  this  con- 
dition occurred.  Normally,  on  subsequent  iterations,  this  point  will 
move  backward  to  the  beginning  of  the  trajectory  and  disappear. 
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V.   PROGRAM  EXAMPLES 

A  single-stage  rocket  trajectory  problem  as  described  below  was 
programmed  utilizing  each  of  the  methods  discussed.  The  ALGOL  computer 
language  was  used  and  the  programs  were  run  on  a  Burroughs  B-5500  digital 
computer . 

Assuming  the  rocket  is  launched  from  an  airless,  non-rotating  Earth, 
the  state  equations  are 

x     =  r  =  V  sin  (y) 
x2  =  6  =  -  cos   (y) 

•    .  *■  _  s  (!-\  /cos(u)   \     »  sin  (y) 

V  §0WUf       "  r2 

\     wo  Isp/, 

£         v  _  g0  /T   \  /     sin(u)    \         \x  cos    (y)        V  cos    (y) 

*     T  wHj*)"  ^     r 

where  r  =  altitude  measured  from  the  center  of  the  Earth,  V  =  velocity, 
y  -   flight  path  angle,  gQ  =  gravitational  acceleration  at  the  Earth's 
surface,  T  =  thrust  (assumed  constant),  u  =  thrust  angle  (measured  from 
the  velocity  vector),  Wq  =  initial  weight,  Isp  =  specific  impulse,  t  =  time, 
[i  =   universal  gravitational  constant,  ©  =  downrange  angle.  The  initial 
conditions  are 

r(0)  =  R  (Earth  radius) 
•8(0)  =  0 
V(0)  =  100  ft/ sec 
7(0)  =  89.87  degrees 

23 


It  is  desired  to  place  the  payload  in  a  circular  orbit  of  maximum  altitude, 
hence 

cp  =  xx(T)  =  r(T) 


*  = 


Appendices  (A)  and  (b)  contain  listings  of  the  steepest-ascent  program 
and  second  variation  program  respectively.   Comments  are  inserted  at  strategic 
points  which  explain  the  sequence  of  operations. 

The  steepest-ascent  program  contains  logic  for  a  single  or  dual  stage 
rocket.   It  was  run  in  the  single  stage  mode  to  generate  a  nominal  control 
history  for  input  to  the  second  variation  program  and  to  compare  results. 
It  was  also  run  in  the  two- stage  mode  to  test  the  coast  optimization  logic. 

The  input  data  for  the  steepest-ascent  program  are 

1.  Initial  velocity  (feet/second)   (must  be  non-zero). 

2.  Launch  angle  (degrees). 

3.  Duration  of  the  first  stage  burn  (seconds),  for  single  stage 
rockets  this  quantity  is  the  total  burn  time* 

h.     First  stage  thrust  (pounds). 

5.  Second  stage  thrust  (pounds),  for  single  stage  rockets  zero 
is  input. 

6.  First  stage  fuel  flow  rate  (pounds/second). 

7«   Second  stage  fuel  flow  rate  (pounds/ second) ,  for  single  stage, 

any  non-zero  number. 
8.  Rocket  liftoff  weight  (pounds). 
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9-   Second  stage  weight  after  separation  (pounds),  for  single 
stage,  any  non-zero  number. 

10.  Initial  value  of  coast  duration  (seconds),  for  single  stage, 
zero. 

11.  Duration  of  the  second  stage  burn  (seconds),  for  single  stage, 
zero. 

12.  Coast  weighting  factor,  A. 

13.  Number  of  stages  (l  or  2). 

The  input  data  for  the  second  variation  program  are: 

1.  Initial  velocity  (feet/ second) . 

2.  Launch  angle  (degrees). 

3.  Duration  of  rocket  burn  (seconds). 

h.     Thrust  divided  by  initial  weight  (pound/pound). 

5.  vx 

6.  v2 

7«   Integration  step  size  and  data  storage  interval. 

8.  K  (See  Section  IV-A) . 

9.  Nominal  control  history. 

Other  parameters  which  the  user  may  desire  to  change  must  be  changed  in- 
side the  program  or  incorporated  into  the  READ  statement . 

For  the  single-stage  runs,  the  following  input  values  were  used: 

Initial  velocity „ . .  .100  ft/ sec 

Launch  angle 89 .87  degrees 

Duration  of  burn 220  seconds 

Thrust 0.0 i+30,000  pounds 

Fuel  flow  rate 1^33.3  pounds/ second 
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Liftoff  weight 333 ,770  pounds 

v1 678,91^ 

V2o -93.58 

Step   size . . . 2   seconds 

K o =  Several  values 

For  the  two-stage  runs,    data  for  the  Gemini-Titan  II   system  were  used: 

Initial  velocity. . .  o 100  feet/ second 

Launch  angle 89.87   to  89. 95   degrees 

First   stage  "burn „ . .  .150   seconds 

First   stage  thrust ^30,000  pounds 

Second  stage  thrust = -  -100,000  pounds 

First    stage  fuel  flow 1666. 6  pounds/ second 

Second  stage  fuel  flow 327  .7  pounds/second 

Liftoff  weight ...331,000  pounds 

Second  stage  weight 70,000  pounds 

Coast   duration .10   seconds 

Second  stage  burn 180   seconds 

Coast  weighting  factor 0.1 

Number  of   stages .2 
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VI.   RESULTS  AND  DISCUSSION 

A.   STEEPEST-ASCENT  VS.  SECOND  VARIATION 

When  the  steepest-ascent  program  was  run  in  the  single-stage  mode, 
the  nominal  trajectory  attained  an  altitude  of  196,015  feet.   The  errors 
in  meeting  the  terminal  constraints  on  the  velocity  and  flight  path  angle 
were  h26   feet  per  second  and  l.lh6   degrees.  On  the  third  iteration  the 
altitude  was  improved  to  260,427  feet  with  terminal  errors  of  O.67  feet 
per  second  and  0.025  degrees.  The  program  was  terminated  when  the  desired 
altitude  change,  dcp,  became  less  than  5,000  f eet  =  At  this  point  fifteen 
Iterations  had  been  completed.  The  terminal  altitude  achieved  was  318,126 
feet  with  terminal  errors  of  0 .6h   feet  per  second  and  0.003  degrees.  The 
associated  control  history  was  punched  on  cards  and  values  of  v  and  v 
were  printed  out.  The  program  ran  five  minutes  and  nine  seconds.   It  is 
estimated  that  this  time  would  be  about  halved  if  the  program  were  run  on 
an  IBM  709O  computer. 

The  output  generated  in  the  steepest-ascent  program  was  used  as  input 
to  the  second  variation  program.  As  expected,  the  nominal  trajectory 
attained  an  altitude  318,126  feet.  On  the  succeeding  forward  integration 
the  trajectory  was  totally  unreasonable.  The  control  deviations  were 
made  smaller  by  increasing  the  initial  magnitude  of  K  but  this  failed  to 
improve  the  results.   Small  variations  were  made  in  the  input  values  of.  v 
which  caused  relatively  large  changes  in  the  results  but  it  was  not  clear 
how  to  make  adjustments  which  would  improve  the  performance  of  the  program. 
The  running  time  was  far  in  excess  of  the  steepest-ascent  program,  taking 
over  three  minutes  to  compile  and  complete  just  one  iteration. 
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B.   TWO  STAGE  ROCKET  TRAJECTORY  WITH  COAST  OPTIMIZATION 

As  stated  in  Section  IV,  the  input  data  for  this  problem  were  those  of 
the  Gemini-Titan  II  launch  system  with  an  arbitrary  choice  of  ten  seconds 
for  the  initial  coast  duration.  Lambda,  the  coast  weighting  factor,  was 
chosen  as  0.1  since  earlier  runs  indicated  that  a  value  of  1.0  caused  the 
coast  variations  to  be  insignificant.  This  choice  proved  to  be  satis- 
factory. 

The  nominal  trajectory  attained  an  altitude  of  392,56V  feet  with 
terminal  errors  of  25  feet  per  second  in  velocity  and  0.84  degrees  in  flight 
path  angle.   On  the  first  iteration,  where  the  program  attempts  only  to 
meet  terminal  constraints,  the  altitude  was  improved  bj'"  18,018  feet,  the 
terminal  errors  were  0.09  feet  per  second  and  0.0025  degrees.  On  the 
fourth  iteration  the  coast  time  was  reduced  to  eight  seconds.  At  this 
point,  a  new  nominal  for  the  portion  of  the  trajectory  following  first 
stage  burnout  was  computed  using  the  new  coast  period.   The  result  of 
this  change  in  coast  was  an  improvement  in  the  velocity  constraint  of 
12  feet  per  second,  and  a  degradation  of  the  flight  path  angle  by  0.25  degrees. 
The  terminal  altitude  achieved  on  this  iteration  was  573,4-50  feet.   On  the 
fifth  iteration  the  coast  period  was  reduced  to  two  seconds,  this  accounted 
for  an  altitude  improvement  of  30,973  feet.  On  the  sixth  iteration  the 
coast  period  was  reduced  to  zero,  this  improved  the  terminal  altitude  by 
13,^77  feet.   In  both  instances  cited  above  where  the  coast  period  was 
changed,  the  terminal  constraint  errors  were  diminished. 

The  altitude  improvement  specified  on  the  second  iteration  was  50,000 
feet.  The  program  was  allowed  to  run  until  this  figure  was  reduced  to  less 
than  1000  feet.  This  proved  to  be  rather  wasteful  as  the  terminal  altitude 
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failed  to  improve  significantly  after  the  desired  improvement  was  reduced 
to  6,250  feet,  which  occurred  on  the  eleventh  iteration.   In  Fig.  1  the 
terminal  altitude,  corrected  for  terminal  constraint  errors,  achieved  on 
each  iteration  is  illustrated.   It  clearly  shows  that  twelve  iterations 
were  sufficient  to  converge  to  the  optimum. 

Figure  2  illustrates  the  convergence  of  the  control  history.  The 
discontinuities  which  occurred  are  due  to  using  discrete  feedback  at 
twenty  second  intervals  rather  than  continuous  feedback.  The  curves  are 
of  different  lengths  due  to  the  changes  in  coast  period  which  occurred. 

Figure  3  illustrates  the  action  of  the  terminal  error  control  scheme. 
Each  time  the  errors  became  excessive  they  were  reduced  to  essentially 
zero  in  one  iteration. 

The  last  terminal  altitude  achieved  was  6l6,573,  a  57*2  per-cent 
improvement  over  the  initial  nominal. 

C.   SENSITIVITY  OF  THE  STEEPEST-ASCENT  PROGRAM  TO  CHANGES  IN  INITIAL 
CONDITIONS 

The  sensitivity  of  the  steepest-ascent  program  was  tested  to  determine 

the  capabilities  of  the  convergence  scheme.  Lauch  angles  of  89=87  and  89*95 

degrees  were  tested  on  the  single  stage  trajectory.  The  resulting  terminal 

altitudes  were  196,015  and  803,428  feeto   In  the  latter  case  the  terminal 

errors  were  1858  feet  per  second  in  velocity,  and  U0.4  degrees  in  flight 

path  angle .   In  both  cases  the  program  converged  in  twelve  iterations  to 

approximately  the  same  optimum  altitude . 

An  attempt  was  made  to  solve  the  problem  given  a  ninety  degree 

launch  angle  which  failed.  Further  testing  at_ lower  launch  angles  was  not 

accomplished  due  to  computer  time  limitations,  however,  it  is  believed  that 
the  tests  conducted  amply  illustrate  the  virtue  of  the  method. 
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VII.   CONCLUSIONS 

k 

The  terminal  control  feedback  scheme  due  to  Bryson  and  Benham  and 

the  method  due  to  Rosenhaum  of  leaving  the  pay-off  unconstrained  in  order 
to  satisfy  terminal  constraints  have  been  combined  and  altered  as  necessary 
to  produce  a  simplified  steepest-ascent  optimization  procedure.  This 
procedure  has  been  shown  to  be  successful  in  solving  a  typical  rocket 
trajectory  problem  including  optimization  of  an  interstage  coast  period. 

In  this  closed-loop  procedure  the  change  in  control  is  computed  at 
each  point  utilizing  continuous  or  discrete  information  on  the  deviation 
from  the  previous  trajectory.   It  is  closed  loop  because  the  program 
continuously  (or  periodically)  checks  on  how  it  is  doing  in  its  attempt 
to  satisfy  the  terminal  constraints.  The  advantage  of  this  procedure  is 
that  larger  deviations  from  the  nominal  can  be  tolerated  while  still 
maintaining  a  bound  on  the  terminal  errors.   It  is,  therefore,  possible 
to  move  rapidly  toward  the  optimum  trajectory  as  is  illustrated  in  Fig.  2. 

The  second  variation  method  due  to  Bullock  has  been  applied  to  the 
same  problem  with  questionable  results.   Consultation  between  the  author 
and  Mr.  Bullock  has  failed  to  uncover  possible  flaws  in  the  theory  or  the 
programming  technique.   Bullock  has  shown  in  several  examples  that  the 
method  is  successful,  however,  in  each  example  the  Hamiltonian  did  not 
depend  explicitly  on  time.   In  the  rocket  problem  this  is  not  the  case. 
Although  this  fact  has  no  theoretical  bearing  on  the  problem,  it  is  the 
only  major  difference  between  Bullock's  examples  and  this  problem. 

Experimentation  revealed  that  the  rocket  problem  is  extremely  sen- 
sitive to  the  choice  of  V.  This  was  not  the  "case  in  the  examples 
presented  by  Bullock.   It  was  thought  that  the  values  of  v  generated  by 
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the  steepest-ascent  program  would  be  very  close  to  correct,  however,  in 
view  of  the  failure,  this  premise  was  laid  open  to  question.  Hence,  one 
reasonable  conclusion  that  can  be  drawn  is  that  V  was  chosen  incorrectly 
and  that  there  is  no  intuitive  or  analytical  method  presently  available 
to  make  the  proper  choice.  There  exists,  of  course,  the  possibility  of 
error,  but  considerable  time  and  painstaking  effort  has  been  expended  to 
minimize  this  possibility. 

Failure  of  the  second- variation  method  notwithstanding,  some  con- 
clusions may  be  drawn  with  respect  to  the  advantages  and  disadvantages 
of  the  two  methods. 

1.  Understanding  the  theory  involved  in  the  second- variation  method 
requires  a  background  in  the  Calculus  of  Variations,  whereas  the  steepest- 
ascent  method  presented  does  not. 

2.  The  second  variation  method  apparently  requires  a  good  estimate 
of  the  sensitivities,  v,  while  the  steepest-ascent  method  only  requires  a 
guess  of  the  nominal  control  and  will  tolerate  a  fairly  poor  guess. 

3.  The  second  variation  method  requires  backward  integration  of  the 
2n2  +  2n  equations  M,  N,  b,  and  A  (where  n  is  the  number  of  state  variables) 
and  a  matrix  inversion  at  every  step  of  the  forward  integration.   The 
steepest-ascent  method  requires  backward  integration  of  $  and  J  which  is 
less  than  2n2  equations  since  J  is  symmetric.  The  matrix  inversion  can 

be  done  at  less  frequent  intervals  using  the  sampled  data  feedback  method 
suggested.   The  second  variation  method  thus  requires  more  computer  time 
per  iteration. 
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VIII.   ADAPTATION  OF  THE  PROGRAMS  TO  OTHER  TYPE  PROBLEMS 

An  attempt  was  made  to  generalize  each  program  so  that  they  could  be 
easily  adapted  to  other  problems.  This  required  a  large  number  of  sub- 
scripted variables  and  matrix  multiplication  loops  in  the  subprograms 
containing  the  differential  equations.   Since  the  subprograms  are  called 
twice  per  integration  step  they  must  be  as  efficient  as  possible.  The  use 
of  subscripted  variables  and  loops  is  most  inefficient  and  results  in  more 
than  doubling  the  running  time. 

Adaptation  to  another  problem  is  still  relatively  simple  however. 
The  programs  contain  sufficient  comments  to  indicate  where  necessary 
changes  must  be  made. 
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APPENDIX  A.   STEEPEST-ASCENT  COMPUTER  PROGRAM 


BEGIN  COMMENT  ERNEST  C.  LUDERS   POX  219   STEEPE S T -  A  SCENT  ; 

INTEGER  FLAG#FA1L#CUUNT,QUITJ 

REAL  T1*C0AST*FF1*FF2*-ISP1*ISP2*WR1>WR2*W1*W2*T2#MM1/MM2# 
*E*  MASQAT/ TOVERWO*  GO,  K,  ISP,  TF*  HH#  SAMPLETIME* 
V0*PREDHFO*PREHFO*DHF0*HF0*J12*J13*PREVF*PREGAMF* 
pKtHF,  DVF*  DGAMF*  30UND*  B0UND2*-  UET,  EX*  CNR*  P* 
MM3/LAMDA*SUM>ULDC0AST>TB2*G>0HFJ 

INTEGER  M A  X I N  OE  X  * L  p M  >  I p  I TE  R  *  $  T  A  GE  S  > OLD  H A  X  I ND  E  X* 
CGMPUTECQAS1*CQASTCOHPUTEJ 

REAL  AKRAY  D  E 1 S  3  * 1 1  3  3  *  A  C  1  r  3  3  #  YP[0l4#0tll0]J 

SAVE  ARRAY  XP/OLDXP C 0 14*0 : 200 J*LP CO : 22* 0  J  200 3*  ERRCO:33* 
T  E  M  P  [  0  *  4  ]  ,    D  E  L  C  C  :  3  3 *  CHKCH3#l«3l* 
TMPC1S2*1?23*  PHI  ,OLCPHI  [0:2003*oEI* 
j  I  N  v/  [  1  :  i  *  1  :  3  3  I 

LA3EL  LI*  L2*  L3*  L4,  L5,L6,HELL*' 


FORMAT  Fl  l//p    X25,  "AT  TIVE  T  =  »p     13*  "   THE  OETERMl", 
"KAN1  OF  J  =  •»* 

LIB. 11*  //*  X33*  3 L2  3.ll*  //*  XI 8*  "CHECK  MATRIX  =  "* 

3  E  2  3 . 1 1 , / / p     v  3  3 ,  3  E 2 3 . 1 1  *  //), 

F2  (X20*  "THIS  RUN  REQUESTED  AN  ADDITIONAL:"*  //)* 

F3  (X2?*F12.2*  "       FEET  OF  ALTITUDE"*  //*  X30*  F9.2, 
"       FEET/SEC  UF  VELOCITY"*  //*  X32,  Ftl.6» 
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"   DEGREES  OF  FLIGHT  PATH  ANGLE**  //)> 
F4  (X20*  "THE  RUN  ACTUALLY  YIELDED  AN  ADDITIONAL*"*  //)* 

F5  (X20*  "TERMINAL  ALTITUDE  =" *  F10.2* 

"        FEET"*  / >     x  2  0  >   "TERMINAL  VELOCITY •SXiZ*"*"* 

XI*  F  9  . 2  »  "        FEET/SEC"*  /  *  X  2  0  * 

"TERMINAL  FLIGHT  PATH  ANGLE         =    H  *  F 1 1 • 6 * 

"    DEGREES"*  /  * 

X  ?  0 »  "ORBITAL  VELOCITY  AT  THIS  ALTITUDE  *  "  *  XI*  F  9  .  2  * 

"        FEET/SEC")* 
F6  CX20*"C0AST  DURA1  IQN",X20*"="*F10 ,2*"        SECONDS")* 
F7  CX20*" CORRECTED  I E  R  M I N A  L  ALTITUDE        ="*F10.2*X7* 

"FEET"), 

HISTURY(4Fia.l2); 

DEFINE  LOGPL  =  FUR  L«-l*2*3*4  DO  t) 

FILE  GUT  CAnQS  0  U.>  10)/ 

COMMENT   MATRIX  MULTIPLICATION  PROCEDURE  GOES  HERE* 


CONVENT   MATRIX  INVERSION  PROCEDURE  GOES  HERE* 


CU^f'ENT-  THIS  PROCEDURE  CONTAINS  THE  SYSTEM  DIFFERENTIAL 
EQUATIONS  AND  IS  USED  On  THE  FIRST  NOMINAL*  ON  RUNS  IN 

which  the  cuast  time  has  oeen  changed  *  and  on  the  final 
precision  run; 
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PROCEDURE  F  U  \  C  T  (  T  >  X,  F > I 

value  t;  real  t;  real  array  x,  ecu; 

BEGIN  WEAL  R>     \l>     $>     C>FEE'TS;  INTEGER  1$    J> 
COMMENT  INTERPOLATE  IN  I  HE  CONTROL  ARRAY  FOR  THE 
PRCPER  VALbt; 

I  «■  tNTIERCT/HH);   J  <-  IF  I<MAXINDEX  THEN  1  +  1  ELSE  II 

FEE  *  CQ  +  PHHI])  ♦  (T/HH  -  I)x(PHI[J]  -  Q)l 

R  «■  X[  l]  +  RE; 

S  «-  S I N  C  X  C  4  3  )  I   C  <-  C  0  S  (  X  [  4  3  )  )   V  «-  X  [  3  ]  I 
TS<-T-T2J  COMMENT  TIME  FROM  STAGE  TWO  IGNITION! 
IF  T  <  T 1  OR  STAGES  =1  THEN  BEGIN 
T  G  V  E  R  *  0  <-  F  F  I  /  w  1 J 

MSRAT*l-TOVErtWOxT/ISPl! 

END  ELSE  TUVERWO<-C; 
IF  T>T?  AND  T2*TF  ThlN  dEGlN 

T0VERwO4»FF2/w2l 

MASRAT4-1-T0VERW0XTS/ISP2I  ENO; 
COMMENT  THESE  AhiE  THE  SYSTEM  DIFFERENTIAL  EQUATIONS/ 

f  c  i  ]  «•  v  x  s ; 

F  [  2  J  <•  vxc/r; 

F  [  3  ]  «-  C Q  *  GOxTOVERWO/MASRAT)xCOSCFEE)  -  K  x  S  /  R  *  2 ; 

f['4j  «■  -k/vxc/r*2  +  qxsin(fee)/v  +  fe2h 
eno  funct; 

comment   this  is  the  backward  integration  procedure  which 
contains  the  adjoint  equations  and  the  controllability 
matrix  equaiions; 


38 


PROCEDURE  LIN8AK  CTd,  LS>  IF))  VALUE  TBJ 

REAL  rb;   REAL  ARRAY  LS,  LFU3J 

BEGIN  REAL  R>  V>     St  Ct  T*  SFt  1.1,  L2t  L3>  L 1 t  IMJ 
INTEGER  ItJ; 
REAL  TS; 

REAL  L5>L6>L7#L8.»Nl>N2>N3,N4#N5>N6>N7>N8j 
REAL  TSM,TCM>M1,M2/M3>FEE*L13>L1A>L33>L34>L43>L44; 
COMMENT  T  IS  BACKWARD  RUNNING  TIME; 

T  «■  TF  -  TB/   I  *  ENTIERCT/HH); 

IF  Iso  then  I  <-  o; 

J    «•    IF     I<MAXINDEX    THEN     1  +  1     ELSE     I;        INT    <■    T/HH    -     IJ 
Ll<-LS[13;       L2<-LSC23t       L  3 «-  L  S  C  3  3 ;       L4«-LSU3; 
L5«-LSt5i;       L6<-LSC6];       L7<-,LS[7]j       L8«-LS[83; 
N1«-LS[9];       N2«-LSU0  3;    N3«-LS[11]I    N4+LSC12  3J 
N5«-L5>C  133;    N6«-LSLH3;    N7«-LS{15  3;    N8«-LSC16  3; 
COMMENT     INTERPOLATE    FUR     THE    PROPER    VALUES    OF    THE     STATES* 
R    <-    C  Q  «■  X  P  C  1 1 1  3  3     «-    I  N  T  x  C  x  P  [  1 ,  J  3    -    Q  )    +    R  E  ; 
V    <-    C  Q  «•  X  P  C  3  p  I  3  3    +    INT*CXPC3>J3-    Q  )  ) 
S    <-     S  I  N  (  (  Q  «-  X  P  t  <i  f  I  ]  )     f     I\Tx(XP[4*J3     -    Q  )  )  ) 

c  «■  co$ccq«-xp[4,i  3)  +  i\tx(xpc4,j]  -  on; 

SF  *  SINC  CQ<-PHI  C 1 3  )  +  INTxCPHI[J3  -  Q)); 
TS<-T-T2t  COMMENT  TIME  TO  GO  UNTIL  STAGE  TWO  IGNITION; 
IF  F>T2  ANU  T2*TF  THEM  BEGIN 
TO.VERr- 0«-FF2/w2; 
MASRAT<-l-TQVERW0xTS/ISP2J  END 
ELSE  TuVERWO<-0; 
IF  T  <T  1  OR  STAGES  =1  THEN  BEGIN 
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T0VERW0*FF1/Wi; 

M ASR AT* 1 "TQ VER WOxT/ I SP 1 J 

E\u; 

COMMENT  THEiE  ARE  THE  AJJUlNT  DIFFERENTIAL  EQUATIONS; 

L F [  1  ]  *-VxCxL2/R*2  *?*KxSxL3/R+3  -  Cx(V  -  2*K/(VxR))xL4/ 

r*2; 
lf r 2 ]  *  o; 

LFC3]  *    +  SxLl  4-CXL2/R  +  Cx(i  +  K/C V*2xR ) )xL4/R  -  GOx 

T0VERw0xSF/CV*2xMASRAT)xL4J 
LFC4]  «■  +Vxcxl1  -  Vx$xL2/H  -KxCxL3/R*2  -  Sx(V  -  K/CVxR)) 
XL4/H; 
LF[5  ]  <■  -VxCxL6/R*2  +2x«xSxL7/R*3  -  C  x  (  V  -  2*K/CVxR))xLa/ 
R*2; 
LFE6]  «•  o; 
LFC?  J  *   t  S  x  L  5  +  C  x  l  6 / R  f  C  x (  l  +  K/(V*2xR))xL3/R  -  G  0  x 

TQV£RW0xSF/CV*2xMASRAT)xL8J 
LF[8  )    <•   +VXCXL5  -  VxsxL6/R  -KxCxL7/R*2  -  5x(V  -  K/CVxR)) 

x  l  a  /  r  ; 

LF[9  ]  «■  -VxCxN2/R*2  +  2xKxSX|\3/R*3  -  C x ( V  -  2xK/(VxR))xN4/ 

r*2; 
LFCIO]  <-  o; 
LFC11]  «■   fSxNl  +C^iJ2/R  +cx(i  +  K/( V*2xR ) )xN4/R  "  GQX 

T0VERW0xSF/CV*2xMASRAT)xN4J 
L  F  C  1 2  ]  «■   +  V  x  C  x  N 1  -  VXSXN2/H  -KxCxN3/R*2  -  S  x  (  v  -  K  /  C  V  x  R  )  ) 

x  N  4  /  R ; 
LF  C  13  3  «■  -VxCxN6/R*2  +2xKxSXN7/R* 3  *  CxCV  -  2XK/C  VxR)  )xN3/ 

r  *  2 ; 

LFC  14  J  «■   0/ 

^0 


LF  C  1  5  J  «■   +  S  x  N  5  +  C  x  N  6  /  R  +  C  x  C  1  +  </(V*2xR))xN8/R  "  GO* 
TQVERW0xSF/CV*2xMASRAT)xN8J 

L  F  [  1 6  ]  «-   +  V x C x N 5  -  V*SxN6/R  -KxCxN7/R*2  *  S  x  (  V  -  K  /  (  V  x  R  ) ) 

<  in  8  /  r  ; 

FEE  *  (Q*PHItIJ)  +  lNIx(PHlC J]-Q); 

TSM  *  GOxTOVE«WO/MASRATxSINCFEE)l 
TCM  *  GOxTOVERWO/MASRATxCOSCFEE)/V^ 
Ml  «•  - L 3 x T S M  +  L a x  r c m ; 
M  2  *  "  N  3  x  T  S  M  t  N  4  x  T  C  M  J 
M  3  *•    -  N  7  x  j  S  M  +  N  8  x  T  C  M  ) 
COMMENT  THEbE  ARE  THE  CONTROLLABILITY  MATRIX  EQUATIONS* 
L  F  [  1  n    *  H  1  *  2  >    L  F  C 1 8  3  «■   M  2  *  2  J  L  F  C  1  9  3 «-   H  3  *  2  J 
LF[20J<-   M1XM2J  LK[21J*   M2XM3!  LF[223*    M3xMlJ 
ENDLINRAK/ 


COMMENT   THIS  IS  THE  FORWARD  INTEGRATION  PROCEDURE  WHICH 
COMPUTES  THL  NEW  CONTROL* 

PROCEDURE  SLOPE  U>     X,    F)i  VALUE  TJ 

REAL  Ti   REAL  ARRAY  X,  F  [  1  3  J 

BEGIN  REAL  R>  V,  S>  C>  T3>  I  NO*  8KIND,  FEE*  TSM,  TCM,  MX, 

\\2>   m j>  t s; 

LABEL  LI.  L2>  L3>  L '4 >     TGOCLOSE * L5, L6 ;     INTEGER  l>    J) 
TS*T-T2;  CUWrtENT  TIME  FROM  STAGE  TWO  IGNITION; 
IF  T  <  T 1  OR  STAGES  =1  THEN  BEGIN 

TOVERWOf-FFl/WU 

MASRAT«-l-TQVERWOXT/ISPi;  ISP«-ISP1J 
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END  ELSE  TOVERW0«-0; 
IF  T>T2  AM)  1211V     THEN  BEGIN 

T0VERW0«.FF2/W2;  ISP*ISP2J 

M  A  S  R  A  T  <■  i  -  T  ij  V  E  ^U)  x  T  b  /  I  S  P  2  ;  E  N  0  / 

IF  T>30llNDx. 999999999  OR  COMPUTECOASTsl  THEN 

BEGIN 
COMMENT  THIS  IS  THE  SAMPLED  DATA  FEEDBACK  SECTION'  AND  IS 
HIT  AT  THE  sampling  internal  and  at  T2  when  a  new   COAST 
period  is  computed; 

if  cqmputecuastm  then 

B  L  (J  I  N 

I  n  o  <-  b  n  u  n  D  /  h  n ; 

H  0  0  w  D  «•  H  t)  U  N  D  +  3  A  M  P  L  E  T  I  M  E  ; 

END  ELSE  IND+ENTIERCT2VHH+. OOOOOOl); 

BKIND«-MAXINDEX-INDJ   SUM<-0  \ 
COMMENT   SINCE  J  GETS  SINGULAR  NEAR  TF*  SKIP  THIS  SECTION 
IF  T  iS  CLOSE  Tj  IF) 

IF  BOUND  >  TF-30   THEN  GO  TO  TOOCLOSE? 

FOR  I«-l>2#3  DO  JEICI'I]  <-  JINV[I*I3«-LPCI  +  16/8KIND]  + 

CIF  T>r2x, 999999  OR  C OMPUTECOA ST= 1  THEN  0[I,I3/ 

lamda  else  o); 

J  lI[  1  ,2  J<- J  IN  VI  l/23«-JEI  [2*  1  3«-JINVC2*  1  3«-LPC20/BKlNDH 

cif  t>t2x, 999999  or  c ompu teco a st  =  1  then  0  c  1  *  2  3/ 
-  la  mo  a  else  o>; 

JET[2/3]<-JiUVL2^3]<-JEir3*23«-JlNV[3,2]<-LPT2t/BKIN0]  + 
CIF  T>T2*. 999999  OR  COMPU TECO A ST= 1  THEN  DC2*33/ 
LAMDA  ELSE  0)> 
JtI[l,33«-JINV[l/3]«-JEIC3*13«-JINVC3*13«-LPC22^BKlN03* 
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(IF  T  >  T  2  *  . 999999  OR  CQMPUTEC0AST=1  THEN  0  [  1  >  3  J  / 

LAMDA  ELSE  0); 

if  CO  -1 P  U  T  E  C  CJ  A  S  T  =  1  THEN  CGMPUTECU  A  ST  +  Q> 

P  <•  T I M  E  C  1  ) ; 

invertcjinv,3M); 

IP  1=1  THEM  BEGIN 

WRITEC<HTHE  J  MATRIX  IS  SINGULAR  AT  TIME  =  H> 
F6.2>>T); 
COMMENT  IF  THE  J  MATRIX  BECOMES  SINGULAR  TERMINATE  THE  RUN; 

gu  ra  hell  end; 

LUQPL  TEMPtL]  *  OLDXPtL* INO 3 $ 

LU.DPL   SUM  «■  SUM  +  LP  [I,     BKINn]x(X[L3  -  TEMPCLDJ 

dllCI]  «■  dhf  -  sum;  sum  <-  o; 

LUOPL   SUM  «■  SUM+LPCl+'8>   BKINOlxCXCLl  -  TEMPCL3)! 

or-H?.)    «■  dvf  -   sum;  sum  «-  o; 

LOOPL   SUM  «■  bUM  +  LpCL+l2>  BKINOlxCXCL]  "  rEMPCLl  )* 
0LLL3]  *■  DGAMF  -  SUM;   SUM  <-  0; 
EKRC1]  *  ERRC2]  «■  ERRC3I  «■  0; 
FUR  L+l'2>3  DO  FUR  M>1,2#3  DO 
ERRCL]  «■  ERRCL]  *  JINV[L*  M]xO£LCMJJ 
CUMMENT   ERK  IS  THE  MATRIX  OF  LAGRANGE  MULTIPLIERS*  THIS 

is  the  end  uf  the  sampled  data  feedback  section; 
tuoclose:  end; 

COMMENT   THE  CONTROL  IS  C  U  1 P  U  T  E  0  BETWEEN  HERE  AND  14 >     NOTE 
THAT  THE  CONTROL  IS  COMPUTED  Ohi.    STEP  AHEAD  OF  THE  CURRENT 
TIME  TO  PERMIT  INTERPOLATION; 
L 1  :    IF  T230UND2X.999999999  THEN 

BEGIN  IF  T=0  THEN  BEGIN  I<-o;   GO  TO  L2  END; 
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L3:       1  <-  B0UND2/HH  +  lj   B0UND2  «■  B0UND2  +  HH; 
IF  80UND2>TF  rHEN  GO  TO  L4J 

LZ:       FEE  «■  PHI [1 3;   V  «■  XPC3>I3;   BKlNQ+MAXINDEX-i; 

TS«.80UND2J 

if   b0und2>t2   1  hen  ts<- dound2-t2 ; 
TSm«-GOxTOVERWO/CQ*(  i-toverwoxts/isp))xsincfee); 

TCH  *  GO*TU\/ErtWO/QxCGSCFEE)/VJ 

Ml  <-  -TSMxLP  L3>3KIND]  +  TCM*LP  [4>BKIND3J 
vd    *    -TSMxLPClWBKINO]  +  TCMxLPT  12*BK  I  NO  ]J 
M 3  «■  -TSMxLPC15*8KINp]  +  TCMXLPU6/BKIND3J 
COMMENT  COMPUTE  THE  NEW  VALUE  CF  THE  CONTROL; 

P  H  I  C  I  ]  «■  (  J  L  «j  P  H  1  [  I  ]  «■  F  H  [  [  I  ]  )  +  M 1  x  e  R  R  C 1 3  +  M  2  x  E  R  R  C  2  ] 

+  M  3  x  E  R  R  [  3  1 J 
If  I=ENTIERCT2/HH+,0000001)+1  AND  STAGES=2  THEN 
CUMPUTECUASTM 
ELSE  COMPUTECOASTfOJ 
IF  T=0  AND  1=0  THEN  GO  TO  L3  END; 
L4!    I  *  LNTIERCT/HH);   J  <■  IF  KMAXINDEX  THEN  1*1  ELSE  II 
FEE  *  C  Q  «■  P  h  I  C  13  D  +  CT/HH  -  I  )  x  (  P  H I  [  J  ]  -  Q  )  ; 

r  <•  a  c  i  ]  +  <  e  ; 

s  «•  siNcxC43);  c  <-  coscx[4]);  v  <-  x[3i; 

CONVENT  HERE  ARE  THE  SYSTEM  DIFFERENTIA!.  EQUATIONS; 

F  [  13  *■  V  x  S  » 

F  C  2  3  <■  V  x  C  /  H  ) 

E[  3]  <-  CQ  «■  GOxTOVERWO/wASRAT)xCQSCFEE)  -  K*S/R*2; 

F[4J  <■  -i$/VxC/R*2  +  9xSlNCFEE)/V  +  FC23; 
COGENT  IN  THIS  SECTION  THE  NEis  COAST  DURATION  IS  COMPUTED; 

IE  CUMPUTECUAST  =  1  UR  CO  A  STCOMP'JTE  =  1  THEN  GO  TO  L5 
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ELSE  GO  TO  L6; 
COASlCOMPUTf>COASTCOMPUTE+lJ 


Aci  3*0; 

A[23*-KxS/R*2  -  F  C  3  3 ; 

A[33«--K/\/xC/R*2  +  F  C  2  3  -  Ft  43i 

8KIND4-MAXINDEX-ENTIERCT2/HH+.  0000001)1 

MM1*A[2]xLPC3*8KIND]  ♦  A[3]x|_P['4>8KINO]J 

MM2*A[2  3xLPCll>3KIrt03  +  A[  3  3  xLPC  1Z#  BK IND  3 J 

MM3*A[23xLPC15>3KIND3  +  A  [  3  1  x|_  P  C  1  6,  8K  I  NO  1 J 

D[l*l]  <•  MM1*25  Dt2>23  <-  MM2*2J  D[3#33  <•  MM3*2J 

D  C  1  *  ^  3  «■  0  C  2  *  1 3  «■  MM1XMK2J 

D[l>3]    <-    D  [  3  *  1  3    «■    f1MlxMM3) 

0  C  2  »  3  3     «■    0  [  3  ,  2  3     *    MM2xMM3> 

IF  C0MPUTEC0AST=1  THEN  GO,  TO  L6? 

CQAS1  <-  COLOCUAS3«-COAST)  +  CQ«-CMMlxERRCl3  +  MM2xERRC2] 

+MM3*£RR[33)/LAMDA); 

COMMENT   CUAST  IS  SET  TU  HE  NEAREST  INTEGER  DIVISIBLE  BY 
THE  STORAGE  INTERVAL*  THIS  AIDS  IN  KEEPING  INDICES  STRAIGHT; 

KRITfc.(<*QELTACOAST  =  n^E20»10>»Q)J 

CQAST«-ENTIER(CC0AST+.9)/2)x2J 

IF  CUAST<0  OR  STAGES  =1  THEN  COASTt-0; 
W  R  I  T  E  (  <  "    COAST  =  fSF2C.lO>>    COAST); 
L6:  END  SLOPE; 

COMMENT-   INTEGRATION  PROCEDURE  GOES  HERE; 


COMMENT  MAIN  PROGRAM  BEGINS  HERE; 

COMMENT  I N 1 1  I A  L I Z  £  bELLb*  FLAGS*  AND  CONSTANTS; 
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ITER^O;   FLAb«-0;   FAIL*  1.5  CQUNT«-OJ  QUIT*OJ 
RE  «■  2.Q9S7;   30  *  32.17;   K  <-  G0xRE*2J 

COMMENT   SAMPLETIME  IS  ME  FEEDBACK  SAMPLING  INTERVAL  AND 
HH  IS  THE  UATA  STORAGE  INTERVAL*   HH  IS  ALSO  THE  INTEGRATION 

step  size; 

sampletime  «•  20  j     hh  *     2   ; 
l 1 !   real'(xp13'0]>xpl4>0j,t1*ff1*ff2'v»r1*wr2>w1#w2>ccast*tb2# 
lam0a*stages)£l43# 

comment   the   input   data   is 

xpt  3>0]   =   ini  i  i  al   velocity 

X  P  [  '1  #  0  ]  =  LAUNCH  A  iv  G  L  E 

Tl  =  FIriST  STAGE  BURN  TIME; 
COMMENT       FFl  =  F1RS1  STAGE  THRUST 

FF2  =  SECOND  STAGE  THRUST 

W R 1  =  FIRST  STAGE  FUEL  FLOW  RATE' 
COMMENT       *R2  =  SECOND  STAGE  FUEL  FLOW  RATE 

M  =  LAUNCH  WEIGHT 

W2  =  SECOND  STAGE  WEIGHT  AFTER  SEPARATION/ 
COMMENT       COAST  =  INITIAL  CHOICE  OF  THE  COAST  DURATION 

T32=  SECOND  STAGE  BURN  TIME 

LAMDA  =  COAST  WEIGHTING  FACTOR; 

comment     stages  =  number  of  stages* 
oldcqasucqast; 

T2*-T  1  +CUAST  ; 

TF«-J2  +  TB2; 

OLOMAXINDEX«-MAXINDEX    <•    ENT  IERCTF/HH+.51  )j 

ISPUFFl/WRi;        ISP2<-FF2/WR2; 

XPC1/0J    <-    XP[2>0]    <-    0;       XP[4,0]    <-    XP[4,0]/57,2957795l3i; 
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COMMENT   GEi\F.RA7E  OK  HEAD  IN  THE  NOMINAL  CONTROL  HISTORY 
FOR  J*G  STEP  1  UNTIL  MAXINDEX  DO  PHI  C 1 1  <■  0* 

COMMENT   INITIAL  CONDITIONS  FOR  BACKWARD  INTEGRATION* 

FOR  J<-2  STEP  1  UNTIL  22  DO  LP[I>C]<-0; 

FOR  1*1  STEP  5  LNTIL  16  0U  LPt  I » 0  3  «•  1 J 
COMMENT   COMPUTE  THE  NOMINAL  TRAJECTORY; 

ADAMSC*'  hH,  0,     IF,     HH#  TF >     Op     C,     XP,  FUNCT); 
COMMENT   P  E  H  F  0  R  M  THE  BACKWARD  INTEGRATION; 
L2:  ADAMS(22,HH,  0,     T  F  j  HH>     TF»  0/  0*  LP  *     L  I  N  9  A  K  )  J 
COMMENT   STurE  THE  OLD  VALUES  CF  THE  STATES  BEFORE  COMPUTING 
A  N E  ft  F 0 R rt A R D  TRAJECTORY; 

FOR  I«-C  STEP  1  UNTIL  MAXINDEX  DO  LOOPL  OLDxP[L>I]  ♦  XP[L>Ii; 
COMMENT   BETWEEN  HERE  ArvQ  L3  THE  ALTITUDE  CHANGE  DUE  TO 
TERMINAL  ERRORS  IS  COMPUTED; 

PREHF  *  XPC  l,MA>INDLXi;   PREVF  <-  XP  C  3>  M  AX  I  NDEX  ]  ; 

P  R  E  G  A  (-1 F  <-  XPE4>MAXINDEX]; 
J12«-LPC20*MAXINDEX]; 
J13«-LPf  22, MAXINDEX  j; 

imp  [  i,  n«-LP  Li  a,  max  i  ndex]; 

rKP[2,13«-TMPn*2J<-LP[21>MAXlNDEXJ; 
TMpt2,23*LPC19, MAXINDEX]; 
INVERTCTMP>2#1>; 

IF  1  =  1  THEN  BEGIN  n K  1  T  E ( <  H  T  H  E  LAMBDA   MATRIX  IS  S  I  N G  U  L  A  R " >  )  ; 

gl  to  iQ   end; 

rnp[i,n*J12xlMP[l#l]+J13xTMP[2MJJ 
rHPCl/2  3«-Jl2xTMP[l,2]+J13xTMPC2*2]; 
VU    «-    SGRTCK/CPREHF+RE)); 
PREnHFO<-Ti,PCl*l]x(VU-PREVF)-TMP[l,23xPREGAMF^ 
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PKEHFD*PRtHF+PREDHFO; 

WKITEC<HOLD  CuRRECTED  HE  =  "  *  F  20  .  1  C>  >  EREHF  0  )  %> 

L3»  IF  FLAG  =0  THEN  BEGIN 
DHF*PREDHF0J 

dvi-  «-sgrtck/cprehf+dhf  +  re>)-prevf; 

DGAMF«--PREfiAMFi  ENDJ 

BOUND  *    B0UNU2  «■  COASTCOMPUTE  *  C  OMPU  TECC  AST  <-  0) 
COMMENT   COMPUTE  THL  NErt  CONTROL  WHILE  INTEGRATING  THE 
SYSTEM  EQUAIIONS  FORWARD; 

ADAMSCA*  Hh,  0*  IF,  HHt     TF,  0*  Op     XP,  SLOPE)/ 
COMMENT   ThL  FOLLOWING  LOGIC  CAUSES  THE  CONTROL  HISTCRY 
TO  BE  PRINTED  OUT  EVERY  FOURTH  ITERATION; 

IF  ENTIER(ITER/«)X4  =  HER  THEN  BEGIN 

WRIT  EC  FOR  1*0  SUP  1  UNTIL  MAX  INDEX  DC  PHICIJ); 

w r i t  e ( [ p  a  g  e  ]  )  end; 
iter«-itlr+i; 

OLDMAXINDEXt-MAXINOEX; 
IF  CQASl^ULDCOASl  THEN 
BEGIN 

COMMENT   ADJUST  STORAGE  LOCATIONS  OF  PHI  AND  FlY  NEW 
NOMINAL  TO  PROVIDE  PROPER  INPUT  OF  STATES  FOR  BACKWARD 
INTEGRATION* 

T2«-ti+coast; 

TF*  T2  +  T<32; 
0L0MAxIN0EX<-MAXINDEX; 

MAXINDEX<-LNTJER(TF/HH+.5n; 
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L«-ENJIER{(CQAST-ULDCQAST)/HH+,0000001); 

FQK    l^LNTIL^fCTl  +  uL^CDAST  )/HH+.00000  0l )     STEP     1 

UNTIL    QlDMAXJNDEx    DO    BEGIN 
PHItl+U  +  PHUlj; 
U     1<E^TIEK(T^/HH+, OOOOOOl)    -1    THEN    PHHI*13*0I 

end; 

I«-£NllER(Tl/HH +.0000001); 

LOOPL       YP£L>03«-XP[L#J.j; 

A D A M S (  4 »    HH>     T 1 1     T  F »     *<-\f    T  F  >       0 »       0  >    Y P t     F U N C T  ) % 

FOR    K*I       STE^     1     UNTIL    MAXINDEX    00 

loopl   xptl>m]<-ypll*m-t.]j 

end; 

COMMENT   PHiNT  OUT  THE  HLSUlTS  OF  THE  ITERATION* 

WRITECE2);   WRITECF3*  DHF>  PVF*   DGAMFx  57, 2957795131 ) J 
WRITECF4);   WRITE(F3j.  XP  [  1,MAX  INDEX  ]  -  PREHF, 

XPC3*MAXIN0Ea]  -  .JRE*/F>  C  XP  1 4, MAXINDEX  3  -  PREGAMF) 

*57. 29577*51 j); 
«RITE  tr5»  XP[1/MAX1NDEX3*  xp[3^aXINDEX]# 

XP[4/MAXlN0EX]x5/'.2957795l31i  SQRT(K/(RE  + 

x  p  [  l  >  m  a  x  i  n  o  e  >:  j  ) ) ) ; 

W R  I  T  E ( F 6 ,  COAST); 

COMMENT       COMPUTE    THE    CORRECTED    TERMINAL    ALTITUDE/ 
VU«-S<3RTCK/(XPC1, MAXINDEX  ]  + RE)); 
UriF0«-TMPll>l]x(V0-XP[3#MAXINDEX])- 

TMPC1/23XXPC4/MAXINDEX3J 
HF0«-XPE1*MAXINDEX]    +    OHFOJ    . 

COMMENT    PRINT    OUT    THE    CORRECTED    TERMINAL    ALTITUDE; 
WKITECF7/HF0); 
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IF  QUIT=J  THEN  GO  TO  16/ 
COMMENT    TLST  TERMINAL  ERRORS  WITHIN  TOLERANCE/ 

IF  ABSlDVF«-CVQ-XP[3/MAXINOEXJ))>50 

0*  ABSCDGAHF«-CXP[4,MAXtNDEX3))  >  .02 
OK  SIGN(DHF)*SIGNCXPC1,MAXINDEX3-PREHF3  THEN 
«EG1N 

FAiL<-FAlL  +  W 
C  u  u  N  T  «-  C  0  U  N  f  ♦  1  / 
COMMENT   TEaT  FOR  IMPROVEMENT  IN  SATISFYING  THE  TERMINAL 
CONSTRAINTS  OR  IMPROVEMENT  IN  THE  CORRECTED  TERMINAL  ALTITUDE/ 
IF  (A8SCDVF)<AtfSCSyRTCK/CPREHF+RE))-PREVF)  ANO 
ABS(DGAMF)<ABS(PREGAMF))  OR  (HFO  >  PREHFO)  THEN 
BEGIN 

KLAG'-i;  DHF<-CIF  'C0UNT  =  1  THEN  0  ELSE  DHFO)/ 
SO  F CI  L5  END  ELSE  BEGIN  FLAG* Of 
FOR  HO  STtR  1  UNTIL  (  i*  AX  I  NDEX«-QLDM  AX  I  NDEX  )  DO 
BEGIN  COMMENT  IF  THE  ITERATION  WAS  UNSUCCESSFUL 

DISCARD  1  HE  DATA  AND  ATTEMPT  ONLE  TO  SATISFY 
TERMINAL  CONSTRAINTS/   ' 
PHi[ I  K-QLQPHIt  IJJ 
LOUPL  XPCL/I  ]<-ULDXP[L/  £]J 
End; 

C  U  A  s  T  <•  0  L  D  C  Q  A  S  T  / 

.  t2«-ti+coast; 

TF*T2+TB2/ 

IF  1TER-1  THEN  BEGIN 
WRITECCPmGEJ  )/ 
COMMENT   IF  THE  F I R S T  ITERATION  IS  UNSUCCESSFUL/  THE 
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FOLIO*  I K1  a  MESSAGE  IS  PRINTED  OUT; 

WRITEC<"THE  TERMINAL  CONSTRAINTS  HAVE  BEEN  VIOM, 
"LATLO  EXCESSIVELY  3EC  AUSE"  >  //>  "THE  NO"* 
"MINAL  CONTROL  OR  THE  INITIAL  CCNOI  T I  0 " # 
"NS  ARE**//, "NOT  CLOSE  ENOUGH  TO  THE  "# 

"OPTimuw guess  again«>);go  rn  L4 

end; 
go  to  l3  end; 

E  N  U 
ELSL  REGIN  FLAG+1*  OriF  <- 1  uOOOOV  (  2*FA  I  L  )  >     COUNTS 
COMMENT   THE  FOLQWING  STATEMENT  CONTROLS  PROGRAM  TERMINATION; 
IF  DHF<10GO  THEN  BEGIN  DHF<-DHFO;  QU I  T<- 1 ;  GO  TO  L5 

end; 
end; 

L  5  »       DVF«-S'JRTCK/(xPLl#^AXlNDEX3*-DHF  +  RE))-XPC3/MAXINDEX]; 

DGAMF«--XP[4>MAXINDEX3; 

GO  TO  L2  ; 
COMMENT   F1'«4L  PRECISION  RUN; 

L  6  I        A  J  4  M  S  C  4  >  i  >     0  #  T  F »     2  0.     20,     $  - 4 ,  a  ~  4  ,  x  P  ,  FUNCT)   ; 
COMMENT  PRINT  QUI  SENSITIVITIES  OF  PAY-OFF  TO  ERRORS  IN 
TERMINAL  CONSTRAINTS; 

w R  it e c  1  m p n , 1 3 • r M p c 1 , 2 ] ) ; 

COMMENT   PUNCH  THE  OPTIMUM  CONTROL  HISTORY  ON  CARDS; 
HELL:WRITECCARDS>HI3TUKY,  FOR  I<-0  STEP  1  UNTIL  MAXINDEX  00 

p  h  i  [  1 3 ) ; 

GO  TO  LU       Lt:   END, 
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APPENDS  B.   SECOND  VARIATION  COMPUTER  PROGRAM 


BEGIN  COMMENT   ERNEST  C.  LUDE&S   BOX  219   SECOND  VARIATION/ 
REAL   RE»GO/K> 1SP#TF , TOV ER WC» NU 1 , NU2> HH, VO*TEl# TE2» 

R*V,TSING#0LDPETM/nETM,0LDTEl»0L0TE2*HUK#T3#Ta# 
OLDNIJ1,  0LDNU2,  OLDT; 
REAL   SUM#FEE>MASRAT>S>CQ,Ci*C2,C3,C4>INT; 
REAL  ARRay  TEMP3»TEMP6Cli33#TEMP4#TEMP5[i:3#ll33# 

DELX* TEMPI*  TEMP2#BCH33#MlNV*M#Ntl!3#lt33# 
C0C0s4#0Jil5]J 
INTEGER   IND,  BKINO; 
INTEGER   ITER>MAXlNDEX#I>ONCE#KK#J>L#STj 
INTEGER  BOUND; 

REAL  ARRAY  OLDXP >  XP , YP C 0  I U >  0 : 1 25 3 , PH I , HU# HUU [ 0 : 1 25 3 > 
FU*HUX[l!3#0»125  3*  7Z[o:24,0:i253> 
FL#SLC1 J  3, 1:3,100 : 1103 *ZPCO: 3,Oi 13; 
LABEL   Ll,L2; 
FORMAT  HlST0RY(ari8.12); 

FORMAT   ri(»ON  THIS  I  TER  AT  I  ON"/  //*  X5,  "N'Jl  =  ", E20.ll* 

"NU2  =  %E20.1l,//,X5#"TEl  =  %E20.H,//> 
X5>"TE2  =  "> 

E20.10,//,X5,"HUK  =  ", E20.ll); 
DEFINE   LOOPl  =  ECR  lM»2>3  DO  U,    LOOPJ  =  FOR  J*-l>2>3  00  t, 
LCOPL  ■  FOR  L«-W2,3  DO  #j 


COMVENT   MATRIX  INVERSION  PROCEDURE  GOES  HERE' 
PROCEDURE  NOMINAL  CT,X/OX); 

VALUE  T;   REAL  If       REAL  ARRAY  X»DxCU; 
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BEGIN 

real  s,c>masrat>fee>g; 

INTEGER  I,j; 

I  *  ENTIERCT/HH);   J  <-  I F  I<MAxlNDEX  THEN  I  +  i  ELSE  IJ 

FEE  «-  (Q«-PHI[I])  +  CT/HH  -  I)X(PHI[J3  -  Q); 

r«-xcij+re; 

MASRAT*1-TOVERW0xT/ISPJ 
S*SINCXU35'   c «-  c  n  5  C  X  C  /4  3  ) ;   v«-xt3]; 
D  X  t  1  ]  «■  V  x  S ; 
DXC2]«-  vxc/r; 
0XC3]  «•  CQ  «■  GqxTOVEpwo/MASRAT)xCOSCFEE)  -  KxS/R*2; 

dx[4]  «-  -k/vxc/"r*2  +  gxsincfeew  +dxc23; 
end  nominal; 

procedure  back(t^z'dz)/ 

value  t;  real  t'  real  array  z>Dzcn; 

BEGIN 

REAL   TB*KS/A#LAMl*LAM2*LAM3/KC*R2*R3>V2/S>C/SF>CF* 

masRat#u#Int; 

REAL    ARRAY    FX/HXX*F#0#SS#M*N#SM»FTN#FM#QNC1:3#H33*CCji 

D*MTD#NTC*HXUC1»3U 
INTEGER    E; 
LABEL    EN; 

TB<-TF-T;       I«-ENT  l£PCT5/hH+. 00 OOOOOOl  d; 
if    i<o    then    i    «-    o;      MASRATi-1-ToVERWOxTB/ISp; 
J    *     IF     KMAXINDEX    THEN    1  +  1    ELSE     IJ        INT    «-TB/HH    -    IJ 


53 


R    «-    CU<-XPU*I3)    +     IN'TxCXPCl,  J]    -    U)     +    RZ> 
V    *    CU*-XPt3*  13  3    +    IN7x(XPC3#J3    -    U)J 

s  *  sin(  clJ«-xp[4,  n )  +   intx(xp[4*j]  -  un; 

C    *■    COSC  (  J«-XPU,  I  ]  3    +     INTx(XPU*J3    -    U))J 
SF    «■    SlNC(U«-PHItl33     +     INTx(PhI[J3    -    U ) )  I 
CF    <-    CCS(  (U«"PHI  [  I  ]  )     +     INTx(PHUJ3     -    11))* 


comment  compute:  partial  derivatives  of   system  equations 
with  respect  to  states; 

F  X  C 1 ,  1  3  «•  0  ;       F  X  CI  >  2  3  «■  S ;       F  X  [  1 ,  3  ]  <-  V  x  c ; 

FxC2,  1  ]«-2xKxS/R*3;       FxC2,2]<-0;       F  X  C  2  ,  3  ]  *--KxC/R*2; 

FX[3^1]«--VxC/R*2  +  2xKxC/(VxR*3); 

FX[3>23*C/R  +  KxC/(VxR)*2-30xTOVe;R^CxSF/CV*2xN'ASRaT)J 

FX[3>3]«--VxS/R  +  KxS/(VxR*2); 

COMVENT   COMPUTE  PARTIAL  DERIVATIVES  OF  SYSTEM  EGU^TlONS 
WITH  RESPECT  To  CONTROL; 

fuci,  n«-o; 

FU[2>I3<--(A<-GOxTOVER^O/MaSRAT)xSF; 

fu[3,  i  ]«-a/vxcf; 


COMMENT  COMPUTE    DERIVATIVES    OF    THE    HAMIlTONIAN; 


LAMUZ[223;       LAM2<-ZC23];       l_AM3«-Z  [243; 
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HUGI3«-LAM2xFU[2#l3+lAM3XFU[3»nJ 

HUU[  I  ]«--i_A^'2xAxCf-LAV3xAxSF/\/    -    HUK; 

COMMENT   CHECK  HUU  >  0/ 

IF  H'JUCI3>0  THEN 
IF  TSI NG=0  TH£N 
BEGIN 

T5I\G«-T8  +  HHx5; 

WRITE(<"HUU    >    0    AT    TIME    =    %F20. 10>#  TB )  I 

WRITeC<"HUU    =    %E20.10>>HUUCI3); 

if  tsing>tf  then 

BEGIN 

tsing*o; 
ts1ng«-1/tsing 

end; 

GO  TO  EN 
END 
ElSE  go  to  en; 

HUX[l,n«-HUX[3,I]«-HXU[l3«-HXU[3]<-0; 

HXUC23«-HUXC2/n«--LAM3xFUC3*I]/v; 

kS«-kxS;   <c«-kxc;  p2«-rxr;  V2«-vxv;  r3«-r2xr; 

HXX[lM]*--6xLAM2xKS/H2*2  +  LAM3x(-6xKC/CR2xR2xV)+2xVxC/R3); 

HXX[i,23*-HXX[2*i]«--LAM3x(2xKC/(V2xR3)     +    C  /  R  2  )  J 

HXX[l,33<-HXXC3/l]«-2xLAf-'2xKC/R3+LAM3x(-2xKS/(R3xV)-(-VxS/R2); 

HXX[2,2]<--LAM3xC-2xAxSF/(V2xv)+2xKC/(R2xV2xV)); 

HXXC2'3]«-HXXC3*2  3«-LAMlxC    +*LAM3x(-KS/CR2xV2)-S/R); 

HXXC3^  33«--LAMlxVxS  +  LAM?xkS/R2  +  LAM3*  (  KC/C  R2xV  )  -VxC/R  )  ; 
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kk«-o; 

loopl  loopj 

BEGIN 

rCL*J3>FXCL*J3-FUCL#I3xHUXCJ*I3/HUUi:n; 

Q[L'J3<-F'J[i,nxFUU>I]/HUU[I]; 

SSCL»  J3<-HXX[L>  J3-HXUCL]       xHXUCJ]       /HUUCI3; 
M[L#J3«-ZCCKK<-KK  +  1)3;       N[L»J]*-Z[KK  +  93; 
COMMENT       STORE    F    AND    SS    MATRICES    FOR    USE    IN    PROCEDURE    NEW|_AMDA; 

if  i >  100  then 

BEGIN 

FICL, J> I3+FCL* J];      • 

S|_CL#J#  I  ]«-SSCL>  J3; 

end; 

end; 

IF    KMAXINDEX    THEN 

DETM«-M[  1,  1  ]x(k'[2,2]xm[3,3]-MC2,3]xMC3,23  )  + 
M[i,23X{M[2,3]xM[3M3-MC2>13xM[3»3J)* 
MC1#33x(mC2*13xmC3*23-MC2*23xmC3#13); 

COMMENT       CHECK    FOR    CHANGE     IN    $lGN    OF    DETERMINANT    OF     MJ 

IF    I<(MAXlNDEX-l )    AND 

S1GmC0LDDETM)/SIGN(DETM)    then 

IF  TSIMG=0  TH£N   BEGIN 

t.sing«-tb+"Hhx5; 

WRITE(<"THE  DETERMINANT  OF  M  CHANGED  SIGN  AT  T  =  "> 
F20,10>#TB);  GO  TO  EN 

END  ELSE  GO  TQ  EN; 
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A«-HU[ I3/HUUCI]; 
LOOPJ  BEGIN 

CCC  J3«--FUC  J»  I3xA; 
Q[  J]«-HXUC  J3  x  a; 

end; 

olcqetm«-detm; 
loopj  loopl  begin 

F  M  i  J  ,  L  3  «■  0  * 

&N[  J/L3*-0' 

sm[  j,l]«-o; 

FTNC  J,L3«"0; 

FOR  KKM*2,3  DO  BEGIN 

FM[J,L3«-FM[J,L3  +  F[J,KK]x  M[KK,L3' 
QNC  J>L3*-QNt  J>  L3  +  Q[J,KK]x  N[KK,L]; 

FTN[J,L]<-FTNCJ^L]+FCKK,J]XNCKK,L]/ 
SM[J,L3«-SM[j,l3  +  SS[J.KK3x  M[KK,L3^ 

end;  end; 
kk<-0;  loo  pl  mtd[l3«-ntc[l3«-qz[l  +  213«-0; 
loopl  loopj  begin 

DZCCKK<-KK  +  l)3<--FM[L>J]+QN[L,Jj; 

DZCkk+93*  SMCL/J3  +  FTNCL>J3* 
DZC(E<-L  +  21)3<-DZCE3  +  FX  C  J,  L  3  *Z  [  J  +  2  1 3  ; 

MTD[L3«-MTDCL3  +  M[J,L]xd[J3; 

ntccl3<-ntccl3  +  nc  j>l  3  xcc  c  j3  * 

end; 

loopl         czcl+1b]«--mtdcl3  +  ntccl3; 
EN:  end  back; 
PROCEDURE  CQNTR0L(T*X>DX)J 
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value  t;  real  t;  real  array  x>Dxcn; 

BEGIN 

LABEL  L1,L2>L3>L4; 

IE  T>ONCE*. 9999999  THEN 
BEGIN 

IE  BOUND=0  THEN 
BEGIN 

lND«-ENTIERC0MCE/HH  +  t0000001); 
GO  TO  L2 

end; 

Li:         lND«-ENTlERCONCE/Hri+.OOOOOOl)  +  i; 
ONCE«-ONCE+HHJ 

L2t         B<lNn<-MAXlNDEX-lND;  IE  I  nD*MAX I NDEX  -  7   THEN  GO  TO  L^j 

kk«-o; 
loopi  begin 

bc i3*zzc i+18»bkind]; 

loo'j  begin 

minv[i*j3<-m[i,j]«.zzcckk«-kk+1)*bkind]; 

MC  I#J]«-ZZtKK  +  9*BKIND] 

end; 
end; 
invertcminv* 3, i )/ 

IE  1  =  1  THEN  WrJTE(<"M  IS  SINGULAR  AT  T  =  ,'>E6.2>>T); 
COMMENT   COMPUTE  INITIAL  CONDITIONS  FOR  PROCEDURE  NEWLAMDA; 
IF  IND=lOO  THEN 
LOOPI  BEGIN 

LOOPJ  BEGIN 
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TEMP4CI,J3«-MINV[I,J3J 

TEMP5tI,J]«-         N[I,J]; 

END/ 

TEMP6Cn<-B[n; 
end; 
if"   ind=102  then 

BEGIN 

delx[13«-xpci>ioo3-oldxp[i»ioo3J 

DELX[2]<-XP[3M00]-0LDXP[3»100]; 
CEL>'C33<-XDt^'100]-nLDXP[a*100]; 

loopi   begin 

sum«-o; 

LOOP  J  5UM«-SUMf.TEMp5CJ#I3xDELXCJ3; 
TEMPI  [  I]«-SUM  +  TEMP6t  I  3; 

end; 

loopi   begin 

sum*o; 

LOOP  J    SUM«-SUM  +  TEMPaCJ,I]xTEMPlCJ3; 

zpc  i>o3«-sum; 

end; 
end; 
comment  compute  coefficients  for  computing  new  control/ 
loopi  begin 
sum«-o; 

loop  j  sum«-suv  +  fuc  j,  ind]xminvc  i,  j3' 
tempi  [  i  ]<-sum; 
end; 
loopi  begin 
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sum«-o; 

LOOPj  SUM«-SU^+TEMPlCJ]xK'[l,j]; 
TEMP2CI  ]«-sum; 

end; 
l^i      if  ind  >maxindex  then  go  to  l3; 

TE^P3Cl]«-0LDXP[l,  inoi; 
TEMP3C2]«-QLDXPC3,  IND]; 
TEMP3t33*QLDXP[ft,IND3; 

T3«-T4«-o; 

LOOPJ  BEGIN 

COCj,IND3*-CHUXtJ*IND3+TE-MP2t  J3  3/HUUCIND3; 
T3«-T3«fTEMPHJ3xZZCJ  +  18*BKIND3; 

end; 

T3«-CT3  +  huC  IND3  )/HUUt  IND  3; 

LOOPJ   T4«-T4  +  C0[  J,  IND3x1EMP3[J]; 

COCa>lNn]«-PHIC  INQ3-TA-T3; 
WR1TEC<"INDEX  =  ",IO,IND); 

WRITEC<"C1  =  %E?0. 10»"C2  =  % E20 .  10 . "C 3  =  M>E20.10> 
»CH    =  M,r:20.10>,  FOR  I«-l  STEP  1  UNTIL  k    DO  CCCIMND3); 
WRITECT3*T4); 
IF  BOUND=0  THEN 
BEGIN 

B0Ur-n«-80UND  +  i; 

GO  TO  LI 

end; 
end; 

L3*     I  *  ENTlERCT/HH);   J  «-  IF  I<MAxlNDEX  THEN  1  +  1  ELSE  11 
IF  TXOLDT  THEN 
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BEGIN 

oldt«-t; 
int*t/hh-ij 

Cl<-(Q«-C0Cl/I]3  +  lNTx(C0Cl*J]-Q); 
C2«-CQ«-C0C2*l3)  +  INTx(C0C2*J]-Q)J 

C3<-CQ«-CD[3,I])  +  lNTx(C0C3#J]-Q); 
C4«-(Q«-C0[4,I])  +  lNTx(C0[4,j]-Q); 

end; 

FEE<-pHnT/HH]<-ClxXm+C2xX[3]+c3xX[a]+CA; 
R«-X[  1  3  +^E'       MASRAT«-l-T0VERW3xT/lSP; 
S*SINCXC43  )'>       C<-CGS(X[/43  );       V«-XC3D; 

o  x  c  i  ]  <-  v  x  s ; 
DXC2J*  vxc/r; 

DXC  3  ]  ♦■  (  3«-GoxTCvERWO/MA5RAT)xCGS(FEE)  -  KxS/R*2' 
0X[a]«--K/Vxc/R*2  +   C*SIN(FEE)/V  +  DX[2]; 
END  CONTROL/ 

COMMENT   INTEGRATION  PROCEDURE  GOES  HERE/ 
PROCEDURE  NEWLAMDA(T>IAM,DLAM); 

value  t;  real  t;  REAL  ARRAY  lam*dlamc l j; 

BEGIN 

REAL  A»INT;Q; 

REAL  ARRAY  D/DELX/SQX   , FTDL C 1 : 3 ] , SDELX [ 1 : 3, 1 00 : 1 10 ] 

>  F  L  I  C  1  J  3  f  1  :  3  ] ; 
INTEGER  IND*KK>Jj; 
LABEL  L1'L2»L3; 
IF  T>ONCEx, 999999999  THEN 
BEGIN 
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Ill 


12: 


1.31 


IF  BOUNO=0  THEN 
BEGIN 

I  ND«-ENTIER(0NCE/HH+. 0000001  ); 

GO  TO  L2 

end; 

ind«-entiercpnce/hh+.0000001)  +  1j 

gnce«-once  +  hh;     if  ind>maxindex  then  go  to  l3; 

delx[1]*-xp[win0]-0ldxpc1mnd]; 

delx[23«-xp[3f  ind3-oldxpc3*ind3# 

delxc33«-xp[4f  ind3-oldxpu'ind3>' 

a«-huc  in03/huuc  ind3; 

LOOpI  BEGIN 

or  I  3«-HUXC  If  INDDxA; 
SDELXCI,  IND3«-0; 

LOOP  J  SDELX[I»lND3*-SDELXLl#IND3+S|_CI»J*IND3x 

DELXC J3; 
SDELXElf  IN03«-SDELXClf  IND3-DCI3; 

end; 

IF  30UnD=0  THEN 
BEGIN 

B0UND«"B0UND  +  1) 

GO  TO  LI 

end; 
end; 

kk«-entierct/hk);  jj<-if  kk<maxindex  then  kk  +  i  else  kk; 
int*t/hh-kk; 
lqopi  begin 

SDXCl3«-CQ«-SDELXCIfKK3)  +  INTx(SDELXrifJJ3-Q); 
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100PJ  FLI[I,J3«-CQ«-FLCI>J>KK3)  +  INTxCn_ri>J,JJ]-Q); 

end; 

l  0  0  p i  begin 

FTDL[I]«-0; 

LOOpJ  FTDLLI  3«-FTDL[  I  J+FLIC  J>  I  3    *LAMCJ]J 

DLAMCI3*   -SOXCll-FTDLCn; 

end; 
end  newlamda; 

comment  main  program  begins  here; 

re«-2.09(?7;   GO<-3?.i7;   k«-g0xre*2;   iter«-o;   lSP«-300; 

OLDT*-6900; 

read(xp[3,0]>xpu,0]>tf,tgverw0>nu1*nu2,hh,huk); 
xp[i>c]<-xfc2'03*-o; 

XPC  n,  Q3«-Xp[4,  03/57.  2957795131; 
MAXINDEX«-EMTIERCTF/HH+.5l); 

READCHISTPRY^FDR  I«-0  STEP  1  UNTIL  ^aXINQEX  DO  PHICI]); 
WRITE(HISTqry#FOR  1*0  STEP  1  UNTIL  MAXINDEX  DO  PHICI3)J 
COMMENT   INTEGRATE  SYSTEM  EQUATIONS  WITH  NOMINAL  CONTROL/ 

ADAMSU*HH,  0*   TF*  HH*  TF>  0,  0,     XP,  NOMINAL); 

V0«-SGRT(K/(X?[1'MAXINDEX3+RE)); 

NU2«-NU2/C2xvO); 

TE1«— XP[A,MAXINDEXI; 

TE2«--V0   +XPC3*MAXINDEX]   ; 
LlJ  IF  ITER/0  AND  ABSC0L.DTE1  XABS(TEl)  AND  A8SC  0LDTE2  }<  A3SC  TE2  ) 

AND  OLDXPC l,MAXlNDEy3>XPC 1 /MAXINDEX3  THEN 

BEGIN 

TEUTEl/2;    TE2<-TE2/2;    HUK«-HUKx  10; 
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end  else  huk«-huk/2; 
comment   compute  initial  conditions  for  backward  integration/ 

r<-xp[i,maxindex3+re;   v«-xpC3#maxinoex]# 

FUR  1*2  STFP  1  OnTIl  18  DO  ZZ[I/03«-0; 

zzc  i>o]«--2xv; 

ZZC4/03  «-   k/r*2; 

Z  Z  [  1  0  /  0  3  +    -4xVxKxNU2/R*3; 

ZZtl2»03  <-  -zzu#<n; 

Z  Z  C  1  3  /  0  3  ♦  -2xKxNU2/R*2J 

Z  Z  [  1  5  /  0  3  *■ + z  Z  C 1  *  0  3 ; 

ZZ  C  17*03  «-  i; 

ZZC  19/  03 «-      -TEi; 

ZZL20/03*    -TE2; 

Z  Z  C  2  1  /  0  3  «■  0 ; 

IF  ITER=0  THEN  BEGIN 

ZZC22'03<-:  +K*NU2/R*2/ 

ZZt23>03*-  2xVxNU2; 

ZZC2a/03«--NUi; 

END  ELSE  BEGIN 

30UND<-0; 

DNCE<-200; 
ADAMSC3,HH,200'TF,TF,TE/^-A,p-4,ZP/NEWLAMDA); 

ZZC22/03«-ZPLl/  13+ZZC22/03; 
ZZ[23/03«-ZPC2/  13+ZZC23»03; 
ZZt26/.0]«-ZPC3/13  +  ZZ[24/0]; 

end; 

tsing*-o;  once*o; 

detm«-o; 
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cogent  integrate  m,n*b*  and  lambda  equations  backward; 

ADAMSC24*  Y*H>     0,     TF>  HH#  10»  0   >0  *     ZZ*8ACK); 

comment  stcre  old  values  of  states  before  running  new  trajectory; 

FOR  I<-0  STEP  1  UNTIL  ^AXINDEX  DO 
FOR  L*l#2*3#4  DO  OLQXPCL*  I  3*XPCL>  1 3  i 

COMSENT   IF  DET(m)  CHANGED  SIGN  THEN  A  CONJUGATE  POINT  EXISTS 
OR  IF  HUU  WENT  POSITIVE  DURING  THE  INTEGRATION  THE  LEGEmDRE 
CONDITION  IS  NOT  SATlSFiED,   TSINg  IS  SET  TO  ThE  TIME  AT 

which  either  occurred  a\d  the  ^orward  integration  is  started 
at  this  point  rather  than  zero/' 

once«-o; 
bound«-o; 

IF  TSING  =  0  THEN  ADAi*S<  ft»  HH,  o»  TF,  HH>  TF,  0>  0*  XP#  CONTROL  > 
ELSE  BEGIN 

IF  TSING>TF  THEN  BEGIN 

WPITEC<"HUU  WENT  POSITIVE  OR  A  CONJUGAL 
"TE  POINT  OCCURRED  TOO  CLOSE  TO  THE  "# 
"END  OF  THE  TRAJECTORYm>);GO  TO  L2  END; 

ONCE<-ENTIERCTSlNG+.l);L«-ENTlERCTSlNG/HH+il); 

FOR  J*-\,2,3,ti    DO  YPtJ*0]«-XP[J#L]; 

ADAMS (4>HH>TSING*TF>HH>HH#0>0,YP> CONTROL); 

FOR  I«-L  STEP  1  UNTIL  MAXINOEX  DO 

FOR  J*l>2>3,4  DO  XP[  J,  I  ]*-YP[  J^I-L]; 

end; 

WRITECFOR  1*0  STEP  1  UnTIl  MAXlNDEX  DO  PHICI]); 

vo«-sqrt(k/cxp£1>  max  index  ]+re))J 
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OLDTEl«-TEi;    0LDTE2«-TE2j 
TEl*-XPt4,MAXlNDEX35 
TE2<--V0   +XP[3*MAXINDEX]   ; 
COMMENT    IF  IMPROVEMENT  IN  TERMINAL  ERRORS  AND  pAYOFF  IS 
LESS  THAN  EPSlLON  (USERS  CHOICE)  THEN  STOP  THE  PROGRAM; 
IF  A35(0L0TE1-TE1  }<.001  AND  A3S ( OLDT E2-TE2  )  <  1 
AND  ABS(OLDXpCl'MAX3NDEX3-XPr WMAxlNDEX]  )<500    THEN 
BEGIN 

WRITECFOR  I*l  STfP  1  UNTIL  MAXlN-DEX  DO  PHICI3)J 

GO  TO  L2 
END  ELSE 
BEGIN 

ITER«-IT£R  +  i; 

WRITE(F1,NU1'NU2,TE3.>  TE2,HU'<); 

GO  TO  LI 

end; 

L25  END. 
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